Defect states and excitations in a Mott insulator with orbital degrees of freedom: 
Mott-Hubbard gap versus optical and transport gaps in doped systems 
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We address the role played by charged defects in doped Mott insulators with active orbital degrees 
of freedom. It is observed that defects feature a rather complex and rich physics, which is well 
captured by a degenerate Hubbard model extended by terms that describe crystal-field splittings 
and orbital-lattice coupling, as well as by terms generated by defects such as the Coulomb potential 
terms that act both on doped holes and on electrons within occupied orbitals at undoped sites. 
We show that the multiplet structure of the excited states generated in such systems by strong 
electron interactions is well described within the unrestricted Hartree-Fock approximation, once 
the symmetry breaking caused by the onset of magnetic and orbital order is taken into account. 
Furthermore, we uncover new spectral features that arise within the Mott-Hubbard gap and in the 
multiplet spectrum at high energies due to the presence of defect states and strong correlations. 
These features refiect the action on electrons/holes of the generalized defect potential that affects 
charge and orbital degrees of freedom, and indirectly also spin ones. The present study elucidates 
the mechanism behind the Coulomb gap appearing in the band of defect states and investigates the 
dependence on the electron-electron interactions and the screening by the orbital polarization field. 
As an illustrative example of our general approach, we present explicit calculations for the model 
describing three t2g orbital fiavors in the perovskite vanadates doped by divalent Sr or Ca ions, such 
as in Lai-a;Sra;V03 and Yi-ajCa^VOs systems. We analyze the orbital densities at vanadium ions 
in the vicinity of defects, and the excited defect states which determine the optical and transport 
gaps in doped systems. 
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I. INTRODUCTION 

Doped vanadium perovskites show a very rich be- 
havior that has been mainly explored by experimen- 
tahsts during the last decade.—"— In contrast to the 
high-Tc cuprates^rJi the perovskite vanadates represent 
quite generally systems with an active orbital degree of 
freedomii^"— The particular excitement with the iJVOs 
perovskite mother compounds, where i?=Lu,Yb,- • • ,La, 
is the quantum interplay of spin-orbital degrees of 
freedom fi2ri2^ which is only under partial control of weak 
Jahn- Teller (JT) and other couplingsi^ The moderate 
coupling to the lattice results from the t2g character of 
the outer valence electrons)^ and explains that the or- 
bital quantum dynamics is not quenched — a property 
vanadates share with other transition metal oxides like 
for example titanates^i^^ and iridatesj-'^--^^ The phase 
diagrams of doped vanadates^ deduced from the resis- 
tivity, specific heat and magnetic measurements show a 
systematic decrease of magnetic interactions and the on- 
set of the metallic behavior only at rather high doping,'^ 
in contrast, for instance, to the high-Tc cuprates where 
the metal-insulator transition (MIT) is found at a few 
percent of hole doping. In Lai-^jSr^^VOs, perhaps the 
most investigated doped vanadium oxide, the MIT is at 
X ~ 0.18 Sr doping)^ while Yi-ajCa^VOa is insulating up 



to X ~ O.SOjA This latter compound undergoes a tran- 
sition from the GC phase to the CG phase^^ at rather 
low doping X 0.02 within the insulating phase. It was 
pointed out^"^ that the double exchange^"— is the mech- 
anism responsible for this transition in antifcrromagnetic 
(AF) phases, similar to the transition from the G-type 
AF (G-AF) to G-type AF (G-AF) order in electron doped 
manganites4i In addition, the change of superexchange 
interaction in the vicinity of the holes^ bound to the de- 
fects has a strong influence on the relative energy of these 
states^S. Remarkably the CG spin-orbital order persists 
up to the MIT. 

The evolution of optical spectra with doping for 
these two vanadium perovskites,'' Lai_a;Sra;V03 and 
Yi-xG&x^Oz, shows that the defects lead to impurity 
states which are responsible for an absorption band deep 
inside the Mott-Hubbard (MH) gap.'*'^ This suggests that 
bound small polarons are the cause of the MIT even at 
high doping concentrations!^ It is eventually the growth 
of the mid-infrared absorption under increasing doping 
and the gradual shift of this absorption toward zero en- 
ergy, which accompanies the MIT. These experimental 
findings call for a theoretical approach that would ad- 
dress the nature of defects and the changes in the mag- 
netic and orbital structure induced by them in a Mott in- 
sulator. Doping in combination with strong correlations 
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leads to spectral weight transfer,- while on the other 
hand the size of the MH gap is not significantly influenced 
by doping.— We discuss these issues in the present paper 
and point out that the defect states here should be distin- 
guished from those observed within the MH gap in other 
circumstances such as the excitonic states in LiCuV04,'*'* 
or the states generated simply by doping an orbitally or- 
dered Mott insulator in absence of defects. 

The description of defect states in systems with 
strongly correlated electrons is an outstanding prob- 
lem in the theory of condensed matter,*^ and has 
been addressed mainly within the Hartree-Fock (HF) 
approach f^^i^"— and the methods based on density func- 
tional theory.—"— It involves for example questions like: 
(i) What is the nature of defect states in a strongly cor- 
related system, i.e., as compared to defects in usual semi- 
conductors or insulators?—"— (ii) What happens to the 
MH gap in the presence of defects? (iii) Which are the 
new features in orbital degenerate MH insulators, as for 
example in the t2g transition metal oxides with active 
orbital degree of freedom, in contrast to doped high-Tc 
superconductors? (iv) Which are practicable methods 
to perform reliable calculations for MH insulators that 
can be extended to take into account defects and disor- 
der? The answer to all these questions is a formidable 
challenge as the defects may strongly affect the subtle 
interplay of several essential degrees of freedom, namely 
spin, orbital, lattice and charge, as we shall see below. 

It is remarkable that the theoretical analysis of dop- 
ing in transition metal oxides, including manganites and 
also superconducting cuprates, is usually performed in 
the metallic regime where the effect of doping on the 
properties of a material with strongly correlated elec- 
trons leads solely to the modification of the electronic 
filling in party filled 3d orbitals.^^ This simplification is 
not possible in the insulating regime where local defect 
fields are not screened, leading to defect states, as for 
example in the case of Ca defects in the lattice of Y 
ions in Yi_2;Ca2;V03,^^ Here, we explore further this 
idea and analyze a model with three t2g orbital fla- 
vors. This is a generalization of the frequently used two- 
flavor model (with active {xz,yz} orbitals) for RVO3 
compounds,^-"— i^^i^Sii^ motivated by earlier electronic 
structure studies?^"— The flrst electron at each V'^+ ion 
occupies the xy orbital, with energy lowered by the crys- 
tal fleld, while the second one occupies one of two de- 
generate orbitals, either the yz or xz, resulting in an 
xy^{yz/zx)^ local configuration. 

The purpose of this paper is to introduce a systematic 
unbiased approach that allows one to study defect states 
in strongly correlated materials. As we show below, the 
interactions, which arise by introducing charged defects 
into a Mott insulator, are numerous and all of them are 
necessary in a realistic approach that is capable of pre- 
dicting experimental behavior of a doped material. We 
adopt here a multiband Hubbard model^ that is sup- 
plemented by several terms responsible for the influence 
of the lattice and defect states on the electronic struc- 



ture. We use it below for the description of t2g elec- 
trons and generalize thereby earlier approaches Jiffii^ 
The electronic structure will be analyzed using the un- 
restricted HF approximation,*'^"''* which is known to be 
reliable in systems with spontaneously broken symmetry, 
for instance for states with magnetic order in multiband 
models for manganites^ Cu02 planes of high tempera- 
ture superconductors^ and low-dimensional cuprates^ 
and may serve to describe complex types of order such as 
stripe or spiral structures 1^"— As we demonstrate below, 
charged defects are the source of impurity states which 
appear as new features within the MH gap, while the 
MH gap itself is essentially unaffected in the low dop- 
ing regime<^ Hence the defect states are responsible for 
the in-gap absorption in the optical conductivity at low 
doping. At the same time, we demonstrate that the HF 
method used here is well designed to treat simultane- 
ously phenomena that arise at distinct energy scales, the 
high energy ~ 1 eV related to the (on-site and intersite) 
Coulomb interaction around the defect states and the low 
energy scale ~ 0.1 eV that is related to the orbital physics 
and controls electronic transport in doped materials. 

The problem addressed here pertains to general tran- 
sition metal oxides with perovskite or layered struc- 
ture, which feature a variety of interesting phenom- 
ena due to strong electron correlations. They in- 
clude MITs in undoped systems, several phenomena 
with long-range (LR) coherence in doped systems, such 
as high-Tc superconductivity,^*'"^'* colossal magnetoresis- 
tance, phases with magnetic and orbital order, and so 
on. In particular, the doped manganese oxides with 
colossal magnetoresistance show several magnetic phase 
transitions^-"^^ that have attracted a broad interest in 
the interplay between spin, orbital and charge degrees of 
freedom in strongly correlated electron systems. In this 
class of systems, orbitals couple strongly to the lattice 
and this coupling supports the superexchange^ lead- 
ing to well separated energy scales for the structural 
and magnetic transitioui^^ In contrast, the i?V03 per- 
ovskites are a challenge for the theory of spin-orbital sys- 
tems as the magnetic and orbital order occur here in the 
same range of temperatures, or even simultaneously as 
in LaVOg^iZi 

A common feature of the i?V03 perovskites is the on- 
set of the G-type alternating orbital (G-AO) order be- 
low the characteristic orbital ordering temperature Too i 
which is in these compounds relatively low. Too — 200 
K, and exhibits non-monotonous dependence on the ionic 
radius of R ionSi-iJ"- The temperature Too comes close 
to the Neel temperature 100 < Tni < 140 K for the mag- 
netic transition to the G-type antiferromagnetic (G-AF) 
phase. This phase competes with another AF phase in 
the i?V03 systems with small radii of R ions.^ In such 
cases, this latter GG phase is stable at low temperature, 
while the former GG phase takes over when tempera- 
ture increases. The best known example of this behavior 
is YVO3 with a transition temperature Tn2 = 77 K. It 
has been shown in Ref. [21I that relativistic spin-orbit 
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interaction is important to describe the reorientation of 
spins associated with the two different types of order at 
this 1st order phase transition. A remarkable reversal 
of the magnetization directio n''^'''^ inside the CG-phase 
with increasing temperature could not be understood in 
the theory until now. 

The principal difficulty in the theory of the i?V03 per- 
ovskites is the quantum interplay of spin-orbital degrees 
of freedom,— which manifests itself in leading contribu- 
tions that emerge due to spin-orbital entanglement.^'*!— 
Quantum effects associated with joint spin-orbital dy- 
namics play a role at zero and at finite temperatures^Zr— 
and the consequences have been observed in several 
experiments:'^* (i) the temperature dependence of the op- 
tical spectral weights,— (ii) the apparent breakdown of 
Goodenough-Kanamori rules for magnetic interactions in 
the CG phase [unlike in other orbitally degenerate sys- 
tems, ferromagnetic (FM) interactions are here stronger 
than the AF ones] ,— (iii) the orbital-Peierls dimerization 
of FM interactions in the C-AF phase22i2ii22 observed in 
the neutron scattering in YVOa,^* and also in LaVOa,^^ 
as well as (iv) the interplay between the spin and or- 
bital correlations, which lead to the phase diagram of the 
i?V03 perovskites.^ The Peierls instability of FM chains 
towards dimerization occurs when spin and orbital de- 
grees of freedom are entangled at finite temperature?^ 
and a previous theoretical study showed that charged de- 
fects favor orbital-Peierls dimerization in the small dop- 
ing regime.^'^ Furthermore orbital fluctuations and their 
competition with orbital-lattice coupling play a crucial 
role for the explanation of the observed non-monotonous 
dependence of the orbital transition temperature on the 
radius of R ions along the BNO^ series. Other 
studies support the picture of strong orbital fluctua- 
tions and their decisive impact on the observed physi- 
cal properties;^"~ We show below that in spite of the 
above quantum nature of the spin-orbital order in the 
i?V03 perovskites, the HF approach is well designed to 
capture the essential features of charge defects. 

The paper is organized as follows. First, in Sec. Ill Al 
we introduce a degenerate Hubbard model for t2g elec- 
trons in the i?V03 perovskites and generalize it to the 
doped materials, such as Yi_3:Caa:V03 with Ca^"*" charge 
defects replacing some Y'^"'" ions. The model includes lo- 
cal and intersite Coulomb interactions, and the Coulomb 
potentials induced by Ca defect states, which change lo- 
cally the electronic structure of t2g electrons within V(3c?) 
orbitals. We describe the treatment of this model in the 
unrestricted HF approximation in Sec. IIIBI In Sec. IH CI 
we present the electronic structure obtained for the un- 
doped CG and GC phases in Mott insulators for the 
typical parameters of RVO3. As explained there, the 
electronic structures provide a realistic description and 
include the multiplet structure of (P excited states. The 
optimized t2g orbitals, at a vanadium ion being the near- 
est neighbor of a Ca defect, change to a new orthogonal 
basis due to the defect-orbital polarization interaction, as 
shown in Sec. IIIII Next, we analyze the electron distri- 



bution on a bond next to a Ca defect in the CG and GC 
phases in Sec. IIV Ai and compare the unrestricted HF 
approach with the exact diagonalization (ED). We also 
analyze one-particle local charge excitations in Sec. IIVBI 
and demonstrate that their energies are well reproduced 
around the defect states when the intratomic Coulomb 
interaction U increases, see Sec. IIV CI Next, we consider 
the electronic structure for the defect states in the CG 
phase in the dilute doping limit in Sec. IV Al and in- 
troduce in Sec. IV Bl the order parameter landscape in a 
doped system. We study the effects of a finite orbital 
polarization in Sec. IV CI and investigate the changes in- 
duced by longer range Coulomb interactions in Sec. IV Dl 
This analysis is followed by the discussion of the elec- 
tronic structure in a correlated material with defect 
states, where we address the gaps observed in the op- 
tical spectroscopy in Sec. IVI Al and the dependence of 
the electron density in xy orbitals on the defect concen- 
tration in Sec. IVI Bl Finally, the paper is summarized in 
Sec. IVIII where we also present an outlook at possible fu- 
ture applications of the general unrestricted HF method 
introduced in the paper in doped transition metal ox- 
ides. The Appendix highlights the importance of Fock 
terms in the present problem and addresses the concept 
of optimized orbitals in the vicinity of a charge defect. 



II. MOTT INSULATOR WITH DEFECT STATES 
A. Degenerate Hubbard model for t2g electrons 

We begin with introducing the multiband Hubbard 
model for t2g electrons designed to describe doped RVO^ 
perovskites, such as the Yi_2;Caa;V03 compounds. This 
effective model includes only t2g orbitals at vanadium 
ions V"^"*", coupled by effective d ~ d hopping elements 
along V-O-V bonds. Below, we consider how this pic- 
ture is modified in presence of charged defects and, in 
particular, when not only the electron density within the 
t2g orbitals of vanadium ions changes, but also the pres- 
ence of Ca defects introduces local interactions acting on 
t2g electrons at vanadium ions in their vicinity. 

The multiband Hamiltonian that describes the MH 
physics of quasi-degenerate t2g electrons and the pertur- 
bations by the defect potentials consists quite generally 
of one-electron terms Hq and two-electron (Coulomb and 
JT) interactions Hint- 

% = Ha +'Hint, (2-1) 
no = Ht + HcF + Hdaf, (2.2) 
Hint = Hu + Hv + HjT. (2.3) 

The one-electron part (|2.2[) is composed of the following 
terms: (i) the kinetic energy {Ht), (ii) crystal-field (CF) 
splitting {Hqf), and (iii) the perturbations generated by 
defects (i?def)- Here, i?dcf includes two terms: 

Hdcf = Hj^p + Hpoi, (2.4) 
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the Coulomb potentials of the charged impurities iJimp 
and the orbital polarization term i?poi- These two terms 
arise due to the presence of Ca defects, and they consti- 
tute effective fields that act on the vanadium ions in the 
neighborhood of defects. They have been introduced in a 
simplified (two-flavor) model with two active t2g orbital 
flavors {yz, zx}^ see below. 

The electron interactions Hint (|2.3p are given by: (i) 
local Coulomb interactions (Hu), (ii) intersite Coulomb 
interactions {Hy), and (iii) the JT effective interactions 
between orbitals which arc induced by lattice distortions 
(Hjt)- In the case of RVO3 perovskites considered here, 
the local Coulomb interactions are quite strong and lead 
to a Mott insulator with high spin (S = 1) and or- 
bital degrees of freedom i^^i^^ Depending on the param- 
eters, the interactions in 'Hint support a particular type 
of symmetry-broken phase?^ and one expects coexisting 
spin and orbital order, either in the form of the CG or GC 
phase, in agreement with experimental observations.^ 

The kinetic energy describes the hopping processes be- 
tween nearest neighbor sites on bonds (ij) oriented 
along one of three cubic directions 7 = a, &, c in the per- 
ovskite lattice. 

Here c]^^ is the electron creation operator at the V'^^ ion 
at site i, with orbital flavor a and spin a —'[, The sum- 
mation runs over all bonds (ij)||7. The effective hopping 
t originates from two subsequent d — p hopping processes 
via the intermediate 0(2p^) orbital along each V-O-V 
bond. It follows from the charge-transfer model with 
p — d hybridization tpd and charge-transfer energy A,'^" 
that t is finite only between two identical t2g orbitals, 
labeled 0(7), that are active along a given bond 
The third orbital, which lies in the plane perpendicular 
to the 7 axis, is inactive as the hopping processes vanish 
here by symmetry. This motivates the convenient nota- 
tion used hereafter^ 

\a) EE \yz), \b) = Ixz), |c) = \xy), (2.6) 

where the orbital I7) inactive along a cubic direction 7 
is labeled by its index. Thus, if we consider an idealized 
case without lattice distortions and without defects, the 
hopping conserves the orbital flavor, 

tlf, = ~U^fi{\~8^o)- (2.7) 

Using the charge-transfer model^S one estimates the hop- 
ping element t — ip^^/A ~ 0.2 eV — this value appears 
to be consistent with that deduced from the electronic 
structure calculations 1^ 

The tig orbital states are non-equivalent in the i?V03 
perovskites due to the GdFeOa-like distortionsj^i and one 
finds that the c orbitals are occupied at every site. This 
is a consequence of the CF splitting term, see Eq. (|2.ip . 



which favors c orbital occupation, 

i 

where = Y.a'^iaa-, with Uiaa = cJacrCjacri IS the elec- 
tron density operator for spin-orbital {acr} at site i. We 
estimated, e° ~ < — hence the c orbital is filled by 
one electron and the second electron at any V'^+ ion oc- 
cupies one of the orbitals in the {a, &} doublet, leading 
to a c^(a, &)^ configuration at each site i. This broken 
symmetry state, which corresponds to electron densities 

(n,c) = l, (n„ = 1, (2.9) 

within tig orbitals at undoped ions, is stable in a 
Mott insulator even for values of Ec as small as 0.5t and 
justifies the two-fiavor model4^ Note, however, that one 
expects that the cubic symmetry with n^-y — 2/3 is re- 
stored at high temperatureSj-iSi^^ but this situation will 
not be analyzed here. 

The defect terms in i?Oi Eq- (|2-2p . describe Coulomb 
interactions of the charged defects and the tig electrons 
within orbitals at neighboring vanadium ions and act on 
them as effective fields. As explained elsewhere when 
an Y^^ ion in YVO3 is replaced by a Ca^^ impurity in 
the doped compound, Yi-^^Ca^^VOs, the lattice is locally 
disturbed and the impurity acts as an effective negative 
charge in the sublattice of Y ions — it is located in the 
center of the Vg cube surrounding the defect site, see 
Fig. 1 of Ref. ,_43l Such a defect is therefore responsi- 
ble for additional Coulomb interaction terms between the 
impurity charge and the electron charges in tig orbitals 
of surrounding vanadium ions. The first additional term 
is the Coulomb potential due to the Ca impurity, iJimp, 
which describes the Coulomb interaction between the de- 
fect at site m with effective negative charge Qm = e (here 
we adopt e = 1) and the total tig electron charge density 
operator at the V ion at site i, 

where Ui = is the electron density operator. 

Here the first sum runs over the Ca defects labeled with 
m € V, where the set V contains all lattice defect sites 
at the considered doping x in Yi-^jCa^^VOa. The sec- 
ond index i labels V ions at distance jr^ — Rm| from the 
considered defect site m. In general, the coordinates of 
the defects {R,„} are statistically distributed, although 
defects will also feel some repulsion that avoids their clus- 
tering. The defect potential of an impurity with negative 
charge Q^i = e is of Coulomb type W{r) = eQo/^cr, 
where ec is the dielectric constant (of core electrons) 4^ 
Below we shall employ a truncated defect potential 

W(|r,-R„,|) = | ^^'if;^^™' (2.11) 
I 0, otherwise 



5 



i.e., i should belong to the cube C^^ = C^(Rm) formed 
by eight nearest-neighbor V sites surrounding the Ca de- 
fect at site TO (Vs cube). With a dielectric constant ec ~ 5 
for YVO3 one finds for the defect potential Vd — 1 eV4^ 
We note that the truncated potential is introduced here 
as it provides lateron a more transparent interpretation 
of the spectra and facilitates the analysis of the energy 
shifts of defect states. 

The last term in Eq. (12. 2p is the orbital polarization 
term, Hpoi- The V ions surround a defect site and form 
a Vg cube^ occupying the positions shown in Fig. [Ija). 
Hpo\ originates from the quadrupolar component of elec- 
trostatic field generated by a charge defect within the 
lattice of V'^^ ions. At each vanadium ion the Coulomb 
repulsion between the defect charge and t2g electrons fa- 
vors the occupation of the linear combinations of t2g or- 
bitals that maximize the average distance of the elec- 
tronic charge from the defect. 

Consequently, the defect-orbital interaction, called be- 
low orbital polarization, takes the form: 



(2.12) 

The coefficients {Aa/3(ri — Rm)} = ±1 depend on the pair 
of considered orbitals a/3 and on the direction — R^; 
they are selected to minimize the Coulomb repulsion. 
This is achieved by 



Aafc (ri — Rm) = 
^ac ij^i Rm) — 
Afcc — Rm) = 



1 li Vi- R„ 
-1 if Vi - R„ 



1 a Vi- R„ 
-1 if r,; - R„ 



1 Vi- R„ 
-1 if Vi - R„ 



(111), (111), 
(ill), (111), 



(2.13) 



(111), (111), 
(111), (111), 



(2.14) 



(111), (111), 
(111), (111). 



(2.15) 



Note that each direction along one of the diagonals of the 
cube involves two vanadium ions. 

The local Coulomb interactions between t2g electrons 
at V^"*" ions are described by;^ 

i.a. i,a<f3 

- 2Jh ^ Sia ■ Sii3 +Jh^ 4at4ai^i0i^il3t- 
i.a.<i(5 i.OL^(3 

(2.16) 

Here the spin operators for orbital a at site i, Sia, are 
defined through the Pauli matrices, {cr^^, cTj^, cr^}, and are 
related to fermion operators in the standard way, i.e., 

^ta = 4a1;'^^»i ' = " n,„;) . (2.17) 



(a) (b) C-AF/G-AO (c) G-AF/C-AO 

(0,1,1) (1,1,1) cTbT ciai ciai cTbT 

(0,0,1 )/T y 



(0,0,0 




ciai I cTbT 
(1,1,0) cM ci 



cToT 




(1.0,0) cib-L cTaT cib-L cTaT 



FIG. 1. (Color online) (a) Reference coordinates (klm) of 
vanadium ions in a representative Vg cube in YVO3. Phases 
with spin-orbital order of perovskite vanadates considered in 
this paper have the electron occupancies in the two magnetic 
phases shown schematically in: (b) (7-AF/G-AO phase; (c) 
G-AF/(7-AO phase. The occupied orbitals {a,6, c} and their 
spin components {t, 4-} are indicated at each corner of the 
cube. 



The first term in Eq. (|2.16l) describes the intraorbital 
Coulomb interaction U between electrons with antipar- 
allel spins; the second one stands for Hund's exchange 
Jh- These two parameters {U,Jh} suffice to describe 
local Coulomb interactions between t2g electrons consid- 
ered here and are determined by the effective Racah pa- 
rameters {A, i?, C};®^ the remaining terms are the in- 
terorbital Coulomb density-density interaction, and the 
so-called "pair-hopping" term which acts in the subspace 
of doubly occupied orbitals. The choice of coefficients in 
Eq. (|2.16p guarantees that the interactions obey the ro- 
tational invariance in the orbital space.-^ Note that this 
Hamiltonian is exact when it describes only one repre- 
sentation of the cubic symmetry group (either t2g or e^) 
for 3d orbitals — this applies to the present case of partly 
occupied t2g orbitals in the i?V03 perovskites. In gen- 
eral, however, the interorbital Coulomb and Hund's ex- 
change interactions are both anisotropioi^ and preserving 
rotational invariance requires also including on-site three- 
orbital and four-orbital terms; more details may be found 
in Refs. IMii 

The LR Coulomb interaction has the usual expression, 



ij ^ia^jP 7 



(2.18) 



where the Coulomb interaction Vij is parametrized by 
the nearest neighbor matrix element Vi 



(2.19) 



with coordinates Vi = {ia,ib,ic) and Vj = (jajbjc), re- 
spectively, given in integer representation. We assume 
that the interaction parameter Vi accounts implicitly for 
the background dielectric function due to the core elec- 
trons. For convenience we define the nearest neighbor 
interaction, 



(2.20) 



in terms of the electron-electron interaction strength Vog, 
where the parameter k = represents the ratio d/a; 
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here a the V-V lattice constant and d is the distance 
between the defect and the nearest V neighbor. In the 
present paper, we discuss results obtained either without 
the intersite Coulomb interaction (V^c — 0) or with the 
full LR interaction. 

The JT term Hjt describes intersite orbital-orbital in- 
teractions that follow from lattice distortions and concern 
the orbital doublet {a, 6}. The distortions are both of 
the JT-type and GdFeOa-type and jointly induce orbital 
alternation in the ah planes and favor identical orbitals 
along the c axis; these interactions are included here by 
the following expression /^'"^^ 

H,T = Vab - E ^^^3 ' (2.21) 

{ij)\\a,b (y>l|c 

where 

T- =]^{n,a-n,b) (2.22) 

is the zih. component of the orbital pseudospin t = 1/2 
along the c axis influenced by lattice distortions. The 
interaction parameters Vab > and Vc > influence the 
subtle balance between the magnetic phases in the i?V03 
perovskites,^ supporting the C-AO order, as realized for 
instance in the undoped YVO3 at low temperatures (in 
the GC phase). These interactions increase within the 
-RVO3 perovskites towards the compounds with small 
ionic radii, and this dependence plays a crucial role in 
the detailed understanding of the phase diagram of these 
perovskites 

B. Hartree-Fock approximation 

The degenerate Hubbard Hamiltonian, introduced in 
the previous Section to describe t2g electrons in a system 
with charged defects, can be solved exactly by ED only 
for very small systems. We present such calculations for 
a single atom and for a bond in Sees. Ill and IV, re- 
spectively. Treatment of larger systems is possible only 
after introducing approximations, either in the strong- 
coupling regime, where charge degrees of freedom can 
be integrated out?iii^ or using HF approximation. We 
adopt here the latter general approach, which requires 
the self-consistent calculation of the full density matrix. 
We performed calculations that describe Yi.^^Caj^VOa 
compound in the regime of low doping using finite clus- 
ters with defects of size N x N x N, and implemented 
periodic boundary conditions (PBC). We present the re- 
sults obtained for 8 in Sees. V and VI. 

The main aim of this paper is to calculate the elec- 
tronic structure of strongly correlated and doped MH 
insulators, i.e., including the defects with their complex 
structure. This implies a computational challenge due 
to the simultaneous appearance of the strong correla- 
tion problem and the perturbations introduced by de- 
fects into the electronic structure. Certainly it requires 



reliable approximations to deal with strong correlations 
and disorder simultaneously. The HF approach is an ef- 
ficient scheme which maps the interacting electron prob- 
lem onto the problem of a single particle moving in a self- 
consistently determined field that stems from all other 
electrons. Whether the HF approach is indeed such a 
reliable scheme that fulfils the requirements stated above 
is the central topic of this paper. 

Before designing this calculation scheme we would like 
to emphasize the following aspects (see also the Ap- 
pendix): 

(i) It has been shown*'^ that a HF factorization when ap- 
plied to a two-orbital model is capable to describe both 
the MH gap and the defect states emerging from the lower 
Hubbard band due to the defects. Even more important, 
the typical multiplet splitting of the transition metal ions 
is obtained if the factorization is performed with respect 
to an optimal local basis set, i.e., with occupation num- 
bers either close to 1 or to 0. That is, the local spin- 
orbital states should be either occupied or empty. 

(ii) This requirement leads to complications in the vicin- 
ity of defects. In particular charge defects lead to a rota- 
tion of the occupied local states, see Sec. IIIIl due to the 
Coulomb potential of a charged defect. 

(iii) In principle, a local rotation can be specified that re- 
moves the off-diagonal terms and defines a new optimal 
local basis, i.e., orbitals having again integer occupancies 
(close to 1 or 0). This simple scheme, however, is not 
successful in presence of kinetic energy terms. Here the 
optimal local or Wannier type orbitals are influenced not 
only by local rotation terms, but by these nonlocal terms 
as well. 

(iv) It is important to recognize that the full unrestricted 
HF scheme, i.e., including all relevant off-diagonal con- 
tributions from interactions, can be used to determine 
the MH-split bands in presence of defects. We will show 
here that this also holds true for the three-flavor case. 
In particular, we will show by comparison with ED that 
the approach provides a surprisingly good description not 
only of occupied states, but also of the unoccupied higher 
multiplet states. 

(v) Moreover, the local rotated basis may serve as an 
intermediate basis and allows one to understand more 
easily the new features emerging due to the defects. For 
example, due to local rotation of t2g orbitals the flavor 
conservation in hopping processes that contribute to the 
kinetic energy is lifted in the rotated basis. That is, in 
the vicinity of defects, orbital flavors mix not only due 
to local orbital rotation (off-diagonal CF terms), but also 
due to the kinetic energy. 

Except for a special attention to all relevant off- 
diagonal Fock terms, the unrestricted HF procedure ap- 
plied here is a standard approach for tight binding mod- 
els of interacting electronsi^"— In general, one decouples 
each term describing a two-particle interaction and re- 
places it by all possible one-particle terms coupled to cer- 
tain effective mean fields — density terms in the Hartree 
approximation (HA), supplemented by off-diagonal Fock 
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terms in the HF approximation. Such expressions are 
augmented by double-counting correction terms, which 
are necessary to avoid double counting of the interaction 
energy. We will adopt just Hartree terms in only two 
cases: (i) when the interaction simulates a potential of 
the lattice and Fock terms would simply be out of scope 
(as for the JT terms), and (ii) when there is no one- 
particle term in the original Hamiltonian that can act as 
a source potential for the related off-diagonal mean fields 
(e.g. with respect to spin degree of freedom). Note 
that the JT terms simulate the potential acting on pairs 
of neighboring ions due to lattice distortions and 
discriminating between the orbital flavors in the doublet 
{a,b}. 

As a matter of fact, Fock terms become active only if 
their off-diagonal mean fields are finite. This can happen 
only if either a one-particle term in the original Hamilto- 
nian (a source for the specific off-diagonal mean field) in- 
duces their finite value, or we choose, as an initial state of 
the self-consistent procedure, a state that requires them 
to be finite. This second case is never realized in the 
studies we perform in this manuscript. Accordingly, we 
will adopt only Fock terms that couple either the same 
orbitals at neighboring sites (as in the kinetic energy) 
or couple different orbitals at the same site (as in the 
orbital-polarization interaction term oc D). These lat- 
ter Fock terms will turn out to be of fundamental im- 
portance to describe faithfully the physics of t2g elec- 
trons in presence of defects. We remark that the HA 
is applied to the on-site interactions Hjj in a way sim- 
ilar to that used to implement local Coulomb interac- 
tions within local density approximation (LDA) in the 
so-called LDA+U method,^ 

Given the above prescriptions, the derivation of the un- 
restricted HF equations is standard and we do not present 
it here in extenso; more details can be found, for instance, 
in Refs. andU^- Following this procedure one arrives 
at an effective one-electron HF Hamiltonian, 



{ij) OLfia 



tafia 

7 t 

C C ■ 
taj fi;a laa jfio 



(2.23) 



This Hamiltonian can be diagonalized numerically, and 
the mean fields appearing in the parameters (see below) 
of the HF Hamiltonian (|2.23l) . and the HF orbitals can 
be determined self-consistently within an iterative proce- 



dure. The orbital energies Siaa are defined as follows, 

£iacr=£c^ca + U (riiaa) + {U - 2Jh) ^(jT-i^ff), 



j{i)eab 



(2.24) 



j{i)<£c 



We have introduced the orbital moments (rj) for the 
nearest neighbors j of the considered site i in the same 
plane ab and along the c axis, respectively, labeled 
The parameter = f if the site i belongs to a cube 
which surrounds a particular defect labeled by m in Eq. 
(|2.1ip . i.e., i £ Cm, and £,im = otherwise. The above 
equations may be further simplified by introducing total 
electron densities {rii} and magnetizations {rrii} at site 
i, and electron densities {riia} and magnetizations {mriia} 
per orbital a at site i, 



m.i = ^{nia^ - n.iai) 

a 

E( 



nii 



(2.25) 

(2.26) 

(2.27) 
(2.28) 



One finds, 



£iaa = £c^ca + - ^Jh^ ""-i - -{U - 5JH)n.ia 

- ^<t{U - JH)rrHa - ^(tJhitIi + ^ Vij{nj) + Vd^^ 



j{i)eab 



(2.29) 



j{i)ec 



Note that the on-site charge repulsive term U{ni — 
^nia) — I Jif (n.j — Tiia) does not contain self-interactions. 
The terms that depend on magnetic moments {m^^} are 
responsible for the magnetic order found in the realis- 
tic regime of parameters, see Sees. V and VI. In the 
weak-coupling regime, where the system is metallic, these 
terms give magnetic instabilities driven by the Stoner 
parameter)2£ being in the present three-orbital model 
/hf = U + 2Jh. 

Near the defect at site m, one finds the local off- 
diagonal elements of the HF Hamiltonian Pia^a at site 
i £ Cm given by the orbital polarization term (I2.12p and 
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by the Fock terms of the Coulomb interaction, 

|3^a^,a^ ^ 6m { i'Ac;. (r^ - R^) - 3 Jff ) (c|^^ c,„„) | , 

(2.30 ) 

where the parameters {^m} are defined as in Eq. (|2.24l) . 
As we show in Sees. Ill and IV, the second term in 
Eq. (I2.30p is crucial for Z? 7^ as it renormalizes the 
orbital mixing term cx D and makes it possible to find 
the orbitals that optimize the energy of the system. The 
hopping parameters Uajiia are renormalized by the Fock 
term that stem from the intersite Coulomb interaction 

(EH, 

tiajfia = ~ ^a^icuy) i (2-31) 

and this renormalization is, for instance, responsible for 
the different bandwidths of the minority and the majority 
bands in transition metals. According to the nearest- 
neighbor nature of the "source" term (the kinetic en- 
ergy) , also the related Fock terms stem from the nearest- 
neighbor Coulomb interaction matrix element Vi. As a 
matter of fact, as discussed above, all other Fock terms 
will simply vanish. 



C. Undoped Mott-Hubbard insulator 

We begin with discussing the reference electronic struc- 
ture of the two phases with broken symmetry in spin and 
orbital space relevant to YVO3: the CG phase and the 
GC phase. The spin and orbital order in these phases is 
shown schematically in Figs. [Hb) and[ljc), respectively. 
Unrestricted HF calculations have been performed for an 
8x8x8 supercell with PBC; larger systems have also been 
considered. In the case of an undoped system, a smaller 
supercell like 6 x 6 x 6 or even 4x4x4 would suffice as 
the HF convergence is indeed very fast and the final re- 
sults for charge, spin and orbital density distributions as 
well as for the total and orbital resolved density of states 
(DOS) do not depend on system size, in both symmetry- 
broken phases, for sufficiently large N [N > 4). Within 
our numerical studies, we have used several parameters 
which are considered realistic for YVO3, but here we 
present only the representative results obtained for two 
parameter sets given in Table ID 

In the relevant regime of parameters for the i?V03 per- 
ovskites, the electrons in undoped YVO3 are localized as 



TABLE I. Standard parameter sets used in the numerical cal- 
culations; all parameters are given in eV. Defect potential Vd 
(|2.1ip contributes only in doped systems. 



set 


t 


e° 


u 


Jh 


Kb 




Vu 


A 


0.0 


0.1 


4.0 


0.6 


0.03 


0.05 


1.0 


B 


0.2 


0.1 


4.0 


0.6 


0.03 


0.05 


1.0 



in a Mott insulator. We consider here the same set of pa- 
rameters as used previously in Refs. [2^ and H^: t = 0.2 
eV, e° = 0.1 eV, [/ = 4 eV and Jh = 0.6 eV. The value 
of t = 0.2 eV has been estimated using the charge trans- 
fer model. For Hund's exchange Jh, we have adopted 
a value that is somewhat screened with respect to its 
atomic value J|f*°™ = 0.64 eV.— Finally, the value of U 
have been selected in order to obtain a value for the su- 
perexchange parameter J = At'^/U = 40 meV, which is 
consistent with the results obtained in neutron scatter- 
ing experiments for spin excitationsi^ As regards the JT 
parameters Vab and Vc, we have chosen values that are 
in the expected range for LaV03. 

The value of the CF potential = 0.1 eV guarantees 
that, in absence of the orbital-polarization interaction 
D, the c orbital is occupied at each V'^'^ ion. In this pa- 
rameter range, and as long as any defect is absent, one 
may simplify the calculation of the electronic structure in 
the HF approximation as the spin of the c electron deter- 
mines the magnetic moment at each vanadium ion, which 
is in the high-spin (HS) state {S = 1) due to finite Hund's 
exchange. In this respect, the HF calculations reduce to 
those performed within the simplified two-flavor model4^ 
Depending on the starting initial conditions, one finds 
two locally stable configurations, the CG and the GC 
phases shown schematically in Figs. [IJa) and[TJb), with 
a uniformly distributed charge of two electrons at each 
vanadium ion. 

The current self-consistent unrestricted HF calcula- 
tions, performed in the undoped regime, give the elec- 
tronic structures of the CG and the GC phase )2i shown 
in Figs. [2] and |3l For the present parameters t2g elec- 
trons are fairly well localized (as t <^ U) and a broken- 
symmetry ground state is found. In such a case, the 
HF approach is rather successful and provides not only a 
good description of the occupied states, but also a quite 
satisfatory description of the excited ones. In fact, in 
the present case also the structure of the upper Hubbard 
band (UHB) is quite well reproduced. To demonstrate 
this, let us consider first the CG phase. At a representa- 
tive atom in position (0,0,0), see Fig. [lib), one finds the 
DOS with well separated occupied states in the LHB and 
with empty states in the UHB, which consists of three 
peaks. The almost-classical ground state corresponds to 
an occupied wave function |c 4, 6 J,). While the c electron 
is perfectly localized {ric — 0.993), the occupancy of the 
b orbital is somewhat lower than one, Uh ~ 0.977. This 
is due to weak fluctuations of the orbital flavor between 
nearest-neighbor sites along the c axis, where all spins 
are aligned (C-AF) and orbitals alternate (G-AO). Con- 
sequently, the I a \.) state has finite low electron density 
{na ~ 0.016). 

The unoccupied part of electronic structure, the UHB, 
consists of three distinct maxima, see Fig. Ufa) . The low- 
est one corresponds to the HS excitation and the system 
is a Mott insulator with a large MH gap: Amh — 2.20 
eV. This value is obtained as the difference between the 
energies measured at the center of gravity (the first mo- 
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FIG. 2. (Color online) Densities of states for the Mott insula- 
tor in the CG phase calculated for an 8 x 8 x 8 periodic cluster: 
(a) total DOS N{u}) (left scale) and the total electron filling 
n(cj) (right scale); (b) local spin-orbital DOSs Niaa{(^) at the 
atom at site i = (0,0,0), orbital a and spin a. The spec- 
tra reflect the spin-orbital order at (0,0,0) and 4,-spin (t-spin) 
states are drawn by solid (dashed) lines, respectively. Vertical 
dashed line indicates the Fermi energy at oj = 0; parameters 
as m set B in Table [B 



FIG. 3. (Color online) Densities of states as obtained for the 
Mott insulator in the GC phase for an 8 x 8 x 8 system: (a) 
total DOS N{u}) (left scale, heavy line), the same for the GG 
phase (thin line) shown for comparison, and the total electron 
filling n{uj) (right scale); (b) local spin-orbital DOSs Niacr{i^) 
at site (0,0,0), orbital a and spin a. Here 4,-spin (t-spin) 
states are shown by solid (dashed) lines, respectively. Fermi 
energy and parameters as in Fig. 2. 



ment) of the peaks corresponding to the highest occupied 
and the lowest unoccupied HF orbitals: \b l) and |a J,), 
respectively, see Fig. ^h). Note that the value of the 
MH gap is given here by the energy of the HS charge 
excitation dfd'j ^ dfdj to the S — 3/2 spin state, 

Eiis = U- 3Jh , (2.32) 

with a a^b^c^ configuration of ^-spin electrons in the fully 
localized (atomic) limit. 

The UHB contains also low-spin (LS) excited states. 
As pointed out before,— , their excitation energies ob- 
tained in the HF approximation (at Ec — 0), 

Slsi ^U-Jh, (2.33) 
Els2 ^U + Jh. (2.34) 

are systematically lower by Jh than the exact ones and 
the upper state is doubly degenerate instead of the lower 



one, as the quantum-fluctuation driven processes (spin- 
flip and "pair-hopping" terms) are neglectedi^ For the 
considered ground state |c | b |), one finds first a LS 
{S — i) state with the \a f) state occupied, see Fig. [DJb). 
The energy of this excitation, measured from the energy 
of the highest occupied \b |) state, is Slsi ~ 3.37 eV, 
which is indeed very close to the value obtained from Eq. 

in the atomic hmit, (f7 - J^) = 3.4 eV. The two 
remaining LS states have energies i?LS2 ~ ^-"^^ ^^'^ 
^LS2 — 4-64 eV, corresponding to double occupancies in 
the c and b orbitals, respectively. These states have been 
labeled with the added electron to the ground state, |c f) 
and \b f) in Fig. [2{b). The above energies are again very 
close to the atomic limit values, ec + {U + Jh) — 4.5 eV 
and {U + Jh) — 4.6 eV, respectively. 

The maxima in the total DOS N{uj) are well separated 
from each other. This can be clearly seen in the total 
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electron filling up to energy uj, 

duj'N{uj') , (2.35) 

-oo 

wliich is almost constant in between the maxima corre- 
sponding to different states, see Fig. [2lja). In particular, 
no states can be found within the MH gap. 

The overall picture obtained for the GC phase is very 
similar, again with the LHB accompanied by three peaks 
in the UHB (see Fig. In this case, the G-AF ground 
state makes the electrons even more localized and one 
finds (down-spin) electron densities ric — 0.993, Ub — 
0.989 and ~ at the representative (0,0,0) atom, see 
Fig. [l](c). The system is a good MH insulator with total 
electron filling rt(0) = 2. A better electron localization 
than in the CG phase may also be concluded from a 
somewhat increased value of the MH gap, Amh — 2.24 
eV. The low-spin excitations are found now at energies: 
£:lsi = 3.57 eV, = 4.53 eV and = 4.62 eV. 

The largest difference with respect to the GG phase of 0.2 
eV is found for the LSI excitation energy i?LSi, see Fig. 
[3Kb). We suggest that this increased excitation energy 
originates from a highly localized (a t) component within 
the ah plane (in the GG phase) as the very same orbital 
is stronger localized here at the neighboring sites along 
a and 6 directions than in the GG phase. We also note 
the increased height of the peaks which indicates again 
stronger localization. 

In summary, we have shown by the HF analysis of both 
GG and GG phases that the DOS consists of the LHB 
and the UHB, and that the latter has a well defined in- 
ternal structure of excited states. The MH gaps are very 
similar in both magnetic structures for the chosen set of 
parameters. In the next sections, we analyze the change 
of local electronic state near charged defects introduced 
by doping and show how such defects change the above 
idealized electronic structure. 

III. OPTIMIZED ATOMIC ORBITALS 



the description of a single atom and compare the HF re- 
sults with those of the ED with the aim at validating 
the scheme that will be used in Sees. V and VI for the 
bulk. The results are mainly presented for varying polar- 
ization parameter Z?, responsible for the adjustment of 
the orbitals to the field which acts on them. 

The set of t2g orbitals {a, 6, c} Eq. (|2.9p . also called 
pristine orbitals hereafter, is selected as orthogonal ba- 
sis in the undoped YVO3 by lattice distortions, which 
lower the c orbitals by a finite CF energy e^. A priori, 
t2g orbitals are expected to adjust to the orbital polar- 
ization field D in the vicinity of a defect and a new set 
of orthogonal orbitals will become the optimal basis for 
the embedded atom. We have determined this optimal 
basis using the self-consistent HF calculations described 
in the previous section and found that all three pristine 
t2g orbitals {a, 6, c} are strongly modified; already for 
Z? = 0.1 eV, the occupied orbitals look quite different 
from {a, 6, c} ones, see below. This gradual orbital rota- 
tion can be easily recognized by analyzing the electron 
densities, n^^ = {n^\), and the off-diagonal elements of 
the density matrix, (cL c^^), shown in Figs. [5] and [51 
respectively, for the GG phase. 

The electron densities are rid = '^H = 1 ^^^^ = 
for the reference state at D = 0. One finds that they 
change fast with increasing D for the orbitals a and h 
(Fig. [S]), while the decrease in the electron density within 
the c orbital is slower. Therefore, in the weak polariza- 
tion regime oi D < 0.02 eV, the two-fiavor model^ is 
sufficient. On the contrary, for D ~ 0.1 eV, all three 
orbitals contribute to the ground state with similar elec- 
tron densities — the electron densities in the a and h 
orbitals are almost equal and ric^ < 0.8. Certainly, the 
two orbitals occupied in this case are quite different from 
the pristine c orbital and one of the doublet {a, b} or- 
bitals and the full orbital space has to be considered. 
This observation is further supported by the values of 
the off-diagonal matrix elements (Fig. [6]), which become 
similar to each other for D ~ 0.1 eV, while for D < 0.02 



To develop systematic understanding of the evolution 
of the electronic structure under increasing doping, it is 
helpful to analyze first the atomic problem for a repre- 
sentative atom near a Ca defect embedded in the GG 
phase of YVO3. Thus, we consider the atomic limit, i.e., 
no hopping (t — 0), as in the parameter set A of Table I. 
The HF mean field terms act on the electronic states of 
such an embedded atom with the ground state |c J, 5 J,), 
see Fig. SJa). This problem is solved self-consistently in 
the HF approximation described in the previous section 
and the solution is systematically compared to the ED 
results obtained for the very same external mean fields. 

Although the electronic configuration cP at an undoped 
V^"*" ion in the vicinity of the defect does not change 
with respect to the Mott insulator described in Sec. Ill CI 
the occupied orbitals are now modified due to the elec- 
tron polarization term Eq. (|2.12l) . Here we focus on 



(a) C-AF/G-AO (b) C-AF/G-AO (c) G-AF/C-AO 
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FIG. 4. (Color online) Reference electronic configurations 
considered in the reference ground state at D = in Sees. 
IIIII and lTVl for: (a) a single atom (dot) in the CG phase; (b) a 
bond (solid line connecting two dots) along the c axis in the 
CG phase; (c) a bond (solid line connecting two dots) along 
the c axis in the GC phase. In each case, the state |c 4- '> 4-) is 
selected as the initial ground state at the (0,0,0) atom. The 
relative position of the defect state in the center of each cube 
is marked by a filled (red) sphere. 
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0.0 0.1 0.2 0.3 0.4 0.5 
D (eV) 

FIG. 5. (Color online) Electron densities n^j, in the pristine 
orbitals 7 G {a, b, c} for an atom embedded in the CG phase, 
see Fig. SJ^a), as functions of increasing orbital polarization 
D. Solid lines from HF, dashed lines from HA and hollow 
circles from ED. Parameters as in set A of Tabled! 




D (eV) 

FIG. 6. (Color online) Off-diagonal elements ^cj^^c^j.^ (7 7^ £,) 
of the density matrix for spin-down pristine orbitals for an 
atom embedded in the CG phase, see Fig. |4ja), as functions of 
increasing orbital polarization D. Solid lines from HF, dashed 
lines from HA and hollow circles from ED. Parameters as in 
set A of Table m 



eV the orbital mixing is large only between the a and 
b orbitals. For the largest value of Z? = 0.5 eV studied 
here, all three t2g orbitals are almost equally occupied 
(Fig. [5]) and the ofF-diagonal orbital elements are almost 
the same for each pair (Fig. [S]). 

It is very remarkable that the HF results presented in 



Figs. [5] and |6] agree exactly with the ED results. This 
shows that the quantum fluctuations, going beyond HF 
approximation, are irrelevant for the ground state in the 
present regime of parameters. In fact, the electron local- 
ization and the symmetry breaking in the spin and orbital 
space, which follow from the large value of [/, suppress 
quantum fluctuations around the exact state with two oc- 
cupied optimal orbitals. In contrast, the results obtained 
within the simpler HA (without Fock orbital terms), turn 
out to be completely unsatisfactory and this approxima- 
tion gives an unrealistic description of the orbital states, 
see Figs. [5] and [51 Within HA, the orbitals b and c are 
almost unchanged and remain still the occupied orbitals 
even for relatively large values of _D > 0.1 eV, where also 
their rotation (mixing) hardly occurs. This demonstrates 
that the Fock terms are essential in the present problem, 
as these off-diagonal terms are responsible for a gradual 
adjustment of the orbital subsystem to the orbital polar- 
ization term Eq. (j2.12p which drives the orbital mixing. 

The gradual evolution of the two occupied rotated or- 
bitals {\b'), |c')} (and the unoccupied orbital {|a')}) rel- 
ative to the original t2g basis as function of the orbital 
polarization field D is illustrated by the overlap functions 
(7 i W i) i'^ Fig. [T] The initial state with occupied or- 
bitals {|&), |c)} changes gradually with increasing D and 
the occupied orbitals {I^OJc')} are linear combinations 
of the pristine ones. One finds once more that, except 
for the region of rather small D < 0.02 eV, all three 
pristine orbitals contribute and the orbitals are strongly 
modified. The new states |c') and \a') arise as linear 
combinations of all three t2g orbitals {a, 6, c}. In con- 
trast, the occupied orbital state \b') has no component 
with |c) orbital character and is just a linear combina- 
tion of {a, b} orbitals. By considering solely this orbital 
state, one could therefore justify a posteriori the two- 
orbital model, "^"^ but we note that: (i) the actual ground 
state involves the occupied |c') orbital too, which has a 
significantly modified shape with respect to the original 
|c) state, and (ii) the presence of finite kinetic energy 
between such orbital states will induce a further redistri- 
bution of the pristine orbital character over the two occu- 
pied orbitals {]&'), |c')}, and will also modify the \a') or- 
bital contributing to excited states. The difference found 
here between the a and b orbital character in the occu- 
pied |c') state follows from the mean fields that arise due 
to the JT interactions. 

The local orbital basis gradually rotates and adjusts 
itself to the orbital polarization when D increases, as 
shown in Fig. [51 The Coulomb repulsion between the 
exceeding electron at the Ca defect in the center of the 
Vs cube, see Fig. [U^a), and the V'^+ ion adjusts the shape 
of the occupied orbitals, {]&'), |c')}, in such a way that 
the electrons involved increase their distance from each 
other, and not too much of the JT and CF energy is lost. 
Therefore, the orbital shapes become all nonequivalent 
for D > and evolve gradually into two geometrically 
similar occupied orbitals {6',c'}, with c' practically lying 
in the plane perpendicular to the direction along which 
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D (eV) 

FIG. 7. (Color online) Fractions of pristine orbitals {I7)} 
in the rotated orbitals {[if')} for atom embedded in the 
CG phase, see Fig. HJa), as functions of increasing orbital 
polarization D. Parameters as in set A of Table HI 
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FIG. 8. (Color online) Optimal orbital basis for an embed- 
ded atom in the CG phase at position (0,0,0), see Fig. |4ja), 
for increasing orbital polarization field D. The defect which 
is the source of the orbital polarization term is located at 
position (i, i, i). The three rotated orbitals are labeled by 
7 — a' ,b' ,c' , with the orbitals b' and c' occupied in the ground 
state. At Z? = eV (first column) {a',6',c'} orbital coincide 
with {a,b,c}. Parameters as in set A of Table [B 



the defect resides, and the empty cigar-like a' orbital, 
directed towards the defect site. This final configuration 
can be easily recognized at the largest studied value D = 
0.5 eV. In the large D regime, the rotated wave functions 
are: 

\a')^±i\a) + \b) + \c)), (3.1) 

\b') = ^{-\a) + \b)), (3.2) 

|c')--^(-|a)-|5)-f2|c)). (3.3) 

In a doped system a hole will go into the topmost occu- 
pied orbital which is an odd linear combination of 
\a) and \b), and does not contain the third c flavor. 

When D is finite and the CF term as well as and the 
mean field terms acting on the atom are neglected, the 
local problem given by the orbital polarization maps onto 
the hole hopping problem on a triangle, with the energy 
spectrum: {—D, —D, +2D}. Having two electrons at the 
V'^"'" ion, the two degenerate states with energy —D are 
then occupied. The increasing similarity of the occu- 
pied states anticipated from this result is indeed observed 
when D increases in presence of other terms. In between 
D = and the largest value D = 0.5 eV studied here, 
one can distinguish two qualitatively different regimes: 
(i) for small D < 0.02 eV only the orbitals {a, b} mix, 
while the c orbital still does not change as it is stabilized 
by the CF energy e° = 0.1 eV; (ii) when D becomes of the 
order of the c orbital is destabilized and all three or- 
bitals rotate further towards their final shapes. The case 
of D ^ 0.05 eV shown in Fig. |5] gives orbitals that are 
already close to those found at the largest studied value 



D = 0.5 eV. Note also that weak JT interactions coun- 
teract the orbital rotation, but the orbital polarization is 
the dominating term for the present set of parameters, 
and the basis rotation is practically completed already 
for D ~ 0.1 eV (see Fig. [T]). 

Orbital polarization infiuences also the one-particle ex- 
citations of an atom close to a defect. Given the ground 
state energy i?o(d") of a Mott insulator in the c?" elec- 
tronic configuration, these excitations in correlated insu- 
lators are calculated as follows: 

— for exciting a hole, as in photoemission (PES) spectra, 

E,n=E„,id''-^)-Eo{r), (3.4) 

— and similarly for adding an electron, corresponding to 
inverse photoemission (IPES) spectra;^ 

£;™ = S„(rf"+i)-So(d"). (3.5) 

We determined the excitation energies for an atom em- 
bedded in the CG phase in proximity of a defect using 
the unrestricted HF approximation and compared them 
with the corresponding ED results, as shown in Fig. |9l 
For the sake of convenience, the chemical potential /i, 
fixed as usual in the center of the MH gap, is subtracted 
in each case. 

In the HF spectra one finds two PES excitations that 
correspond to removing an electron either from b' or c' 
orbitals; the energy difference between them is small and 
comes from the CF splitting £°. These states form the 
occupied LHB. The apparent decrease in energy present 
in Fig. [51 {dE„i/ dD) = —1.5 instead of the expected 
{dEm/dD) = —1, follows from the subtraction of the 
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FIG. 9. (Color online) Excitation energies Em for an atom 
embedded in the CG phase, see Fig. HJa), measured with 
respect to the chemical potential, as functions of increasing 
orbital polarization D. Solid lines from HF (heavy/light lines 
for occupied/unoccupied states); and hollow circles from ED 
(the color scale on the right indicates the spectral weight of 
each level, as well as the size of the symbols). The MH gap 
is found between the LHB and the HS state of the UHB; the 
labels of the LS excited states refer only to the regime D ~ 0. 
The chemical potential /i is fixed at the center of the MH gap 
Amh. Parameters as in set A of Table [H 
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FIG. 10. (Color online) Excitation energies Em for an atom 
embedded in the CG phase, see Fig. |4ja), obtained in the 
HF (solid lines) and in the HA (dashed lines) for increasing 
orbital polarization D. The energies Em are measured with 
respect to the chemical potential /i (dotted line). Parameters 
as in set A in Table H] 



chemical potential. The four excitations at positive en- 



ergies correspond to the four possible (P excited states 
analyzed for the undoped case in Sec. Ill CI (i) the HS 
state has the lowest energy and the \a 1) state is occu- 
pied, (ii-iv) in the three LS states with higher energies a 
f-spin state is occupied in one of the three orbitals. No- 
tably, these excitation energies depend rather strongly 
on the orbital polarization D. The energy of the HS 
state [S^ = — |) increases with increasing D when all 
three orbital flavors are occupied with the same spin, 
and (dE^/dD) = 2 before chemical potential subtrac- 
tion. As a matter of fact, such increase concerns also the 
LS2 excited state with = and |c' f) state occu- 
pied. On the contrary, the remaining two LS states of 
the UHB have either (a' f) or (6' t) occupations by the 
t-spin and their energies decrease. Note that the energies 
of the HS and of these latter LS excited states cross for a 
value oi D > 0.5 eV, but we have estimated that so large 
values of D are unrealistic. 

A comparison of the HF results with the exact results 
found using ED for an embedded atom close to a defect 
gives an excellent agreement for the LHB and even for 
the lowest HS excited state of the UHB. This may be ex- 
pected as these three states have no quantum corrections. 
On the contrary, as already pointed out in Sec. IIIB| in 
the LS sector one finds that the excitation energies for 
D = are lower by ^ Jh than the exact values 4^ When 
D increases, the spectral weight moves from the high- 
energy LS2 state with — /i > 4 eV to the lower energy 
excitations with Em — fJ. ^ S eV. This is in contrast to 
the HF results, where the one-electron wave functions do 
not allow for any spectral weight transfer, and instead 
the excitation energies change and interpolate between 
the relevant energies in a continuous way. Already for 
D ~ 0.15 eV, one finds that the spectral weight found in 
the ED has moved completely to the lower energy regime 
characteristic of LSI states, and when D increases fur- 
ther there is one excitation with large spectral weight 
that follows roughly the two lower HF excitations of LS 
character, and the third one with increasing energy, sim- 
ilar to the one found in the HF, see above. Rather small 
energy difference between the ED energies and the HF 
ones found here demonstrates that the quantum effects 
are quenched in the orbital system when the orbital po- 
larization D is sufficiently large. Even more importantly, 
a satisfactory agreement between the energies found in 
the HF approximation and the exact values in the LS 
sector shows that the HF states correctly adjust them- 
selves to the underlying interactions and simulate the 
actual multi-electron states in a very realistic way. 

One may wonder again whether the HA would not give 
at least a satisfactory description of the excitations being 
close to the MH gap, as we have shown above that for 
all of them the quantum corrections to the HF energies 
are negligible. Instead, the HA fails also here in a rather 
spectacular way — one finds that the MH gap exhibits 
a qualitatively incorrect behavior and decreases with in- 
creasing orbital polarization D, and even the structure 
of the LHB, with two PES excitations of almost equal 
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energies, is not correctly reproduced, see Fig. [10] Large 
systematic deviations are also found for the higher LS 
excited states. All these misleading features follow from 
the incorrect local basis rotation with increasing orbital 
polarization Z?, and thus from the incorrect evolution of 
the orbital electron densities shown in Fig. \5\ It can 
be concluded that, at least in case of conflicting external 
fields (e.g. the CF, the JT interaction and the orbital 
polarization), the HA is totally unsatosfactory. 



IV. A BOND NEAR THE DEFECT STATE 

The electronic structure in the vicinity of defects de- 
pends also on the kinetic energy (|2.5p . While in the un- 
doped bulk (Sec. IIIC|) . the kinetic energy is diagonal and 
couples two atomic t2g states with the same orbital flavor 
at neighboring sites, this is not longer the case close to 
defects where the local orbital basis has to adjust to the 
actual fields acting on each atom, see Sec. IIIII There- 
fore, one expects that qualitative changes may arise with 
respect to the case of an embedded atom when electron 
hopping contributes. The analysis of an embedded bond 
presented in this section addresses this issue and serves 
to deepen and improve the understanding of the HF re- 
sults presented for the bulk in the low doping regime, 
analyzed in Sees. V and VI. 



A. Electron densities induced by orbital 
polarization 

We consider first a representative bond {ij) in the CG 
phase, shown in Fig. SJ^b). The c orbitals are occupied 
at both sites of the bond (in the absence of defects, i.e., 
at I? = 0) and the JT interactions (|2.2ip favor the b (a) 
orbital occupation at the lower (upper) site. This case is 
the most interesting one as the spin order is FM and the 
two electrons in the {a, b} orbital doublet can delocal- 
ize along the bond. The kinetic energy is strongly sup- 
pressed by large local Coulomb interactions, but still one 
finds, within the ED, that orbital fluctuations along the 
bond contribute and modify the electron densities riij^ 
for 7 = a, & at = 0, see Fig. [TTJ Such fluctuations are 
very important and support the FM coupling along the 
vertical bonds in the undoped material.>i2. As we show be- 
low, the HF approximation may be still used although it 
is unable to treat orbitally entangled states^^i^^ as these 
fluctuations are quenched by external flelds, including 
orbital polarization D > 0. In contrast, no fluctuations 
occur for the localized c electrons and their density is 
perfectly reproduced by the HF approximation. 

As in the case of a single atom analyzed in Sec. IIII[ an 
increasing orbital polarization D modifies the occupied 
orbitals and thus the electron densities, see Fig. [TTJ Full 
symmetry in the orbital occupancies is preserved, pro- 
vided one interchanges the orbital flavors {a, b} -H- {6, a} 
between the sites i = 1 and j = 2. Similarly as in the 
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FIG. 11. (Color online) Occupations of pristine orbitals n^^j, 
at the embedded bond (ij) along the c axis in the CG phase, 
see Fig. |3Ib), as functions of increasing orbital polarization D. 
Solid lines from HF and hollow circles from ED. Parameters 
as m set B in Table HI 



embedded-atom case, the orbitals of the doublet active 
along the c axis mix easier and this mixing competes 
with the JT effective interactions, which act on them at 
each site. The orbital polarization overrules the effective 
molecular fields, which act on the {a, b} orbitals, already 
for D = 0.03 eV (see Fig. [IJ), and the trend in the ED 
results changes — in the regime of small D, the filling of 
the occupied {1&, 2a} orbitals, nn,^ — n2a^, first slightly 
increases with increasing D with respect to the one at 
-D = 0, but next when D increases further it starts to 
decrease and the orbital entanglement is suppressed. In 
this latter regime, one again finds a perfect agreement 
between the HF and the ED results, which shows once 
more that the HF method is very reliable for states with 
broken symmetry. Note that, as long as the orbital fluc- 
tuations dominate for very small values oi D < 0.02 eV, 
the orbital electron densities almost do not change, being 
stabilized by the orbital entangled state that also coun- 
teracts the rotation of c orbitals in this regime. On the 
contrary, in the regime of large D > 0.1 eV, the electron 
densities are very similar to those found in Sec. IIIII for 
an embedded atom, cf. Fig. [5] 

It is of interest now to compare the evolution of elec- 
tron densities with increasing D between the CG and the 
GG phases. The latter phase features much less quan- 
tum fluctuations, as the AF spin order is present here also 
along the c axis and flnite electron hopping would gen- 
erate LS states unfavored by Hund's exchange. There- 
fore, the hopping is suppressed also along vertical bonds 
in agreement with the double exchange mechanism, at 
least in the HF picture where quantum fluctuations are 
almost completely neglected. Therefore, even when the 
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FIG. 12. (Color online) Occupations of pristine orbitals n^^j, 
at the embedded bond (ij) along the c axis in the GC phase, 
see Fig. Oc), as functions of increasing orbital polarization D. 
Solid lines from HF and hollow circles from ED. Parameters 
as m set B in Table HI 
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FIG. 13. (Color online) Excitation energies Em for the embed- 
ded bond (ij) along the c axis within the CG phase, see Fig. 
UJb), as functions of increasing orbital polarization D. Solid 
lines from HF (heavy /light lines for occupied/unoccupied 
states), and hollow circles from ED (the intensity scale and 
marker size gives the spectral weight of each level). Parame- 
ters as in set B in Table HI 



electron hopping is finite, the kinetic energy is totally 
suppressed in this symmetry-broken state at Z? = 0, and 
the electron densities are as in Eqs. ()2.9p . see Fig. [T2] 
When D increases, the orbital rotation towards the op- 
timal orbital basis in the limit of large D takes place, 
see Fig. [51 and the orbitals are again locally adjusted by 
the orbital polarization term. Notably, the electron den- 
sities for the initially filled {16,2a} and empty {la, 26} 
orbitals are here more distinct than in the CG phase 
as the kinetic energy, which helps to make the electronic 
distribution almost symmetric, is absent. Altogether, the 
excellent agreement found in this case between the HF 
and the ED results confirms the more classical character 
of the GC phase with respect to the CG one. 



B. One-particle excitations 

We determined also the excitation energies, defined in 
Eqs. p.4p and p.Sp . for the embedded bond. Overall, 
the results found for the CG phase, see Fig. [131 resem- 
ble those obtained for an embedded atom, see Fig. [HI 
A qualitatively new feature in the ED is the broaden- 
ing of the Hubbard subbands, found for the LHB and 
for the individual excitations that belong to the UHB. 
This broadening comes from the action of the kinetic en- 
ergy which mixes the occupied states at both sites of 
the bond (ij) leading to bonding and antibonding states 
that are modified by the JT effective interactions. As 
for the atom, U is the dominating energy scale aand the 
MH gap increases linearly with D. There are four excita- 
tion energies that belong to the LHB and stem from two 



excitations at each ion of the considered bond. Having 
these splittings in the one-electron states, it is natural 
to expect that also the three-electron states will occur in 
certain energy intervals for each (HS or LS) excited state. 
The number of excited states in the HS states is two as 
there could be two states filled by three electrons, one at 
each site. 

The analysis of the LS states is more subtle here than 
in the embedded atom case, as not only four LS2 states 
(two at each atom) split off with increasing D, but be- 
cause also their spectral weight is transferred from the 
LS2 states to energies characteristic of the LSI excita- 
tions. All in all, the qualitative picture found before for 
the embedded atom is confirmed here for the embedded 
bond, and again the LSI excitations, which correspond to 
creating doubly occupied orbitals of 6 character or single 
occupations in the a f orbital, stay roughly in the same 
distance from the LHB independently of the value of D, 
while the LS excitations at higher energies are finally 
the ones that involve double occupancies in the c orbital. 
The distance between these states and the HS excitations 
is close to 3Jh, the value deduced from the multiplet 
splitting.^ We remark that the HF approximation gives 
excitation energies that reproduce quite well the values 
obtained within the ED. This holds in particular for the 
LHB and for the HS states of the UHB. 

Given the more classical character of the GC phase, 
the LHB and the UHB are even more pronounced in this 
case, see Fig. [Ml Here, the splittings arising from the fi- 
nite kinetic energy and found in the CG phase are absent. 
Accordingly, the dependence of the one-hole (PES) and 
the one-electron (IPES) excitation energies is very similar 
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FIG. 14. (Color online) Excitation energies E,n for the embed- 
ded bond {ij) along the c axis within the GC phase, see Fig. 
m^c), as functions of increasing orbital polarization D. Solid 
lines from HF (heavy/light lines for occupied/unoccupied 
states), and hollow circles from ED (the intensity scale and 
marker size indicate the spectral weight of each level). Pa- 
rameters as in set B in Table |T] 



to those found for a single embedded atom in this phase 
(not shown). The LHB consists of two excitations that 
correspond to adding a hole to one of the occupied or- 
bitals — either the c orbital, at the lower energy because 
of the crystal field, or the orbital in the {a,b} doublet 
locally favored by the JT interactions. These excitations 
have been found both in the ED and in the HF at exactly 
the same energies. 

The situation for the electron excitations at high en- 
ergy is again more subtle. Here one finds in ED four 
LS excitations at energy U with respect to the ground 
state, and they split off with increasing orbital polar- 
ization D — the ones with a roughly constant energy 
distance from the LHB do correspond to the double oc- 
cupancy of b orbital or single occupations in the a f or- 
bital, while the ones that follow from double occupancy 
of c orbital exhibit a similar energy increase with increas- 
ing D as the HS state. Already in the range of relatively 
small D ^ 0.05 eV, the spectral weight is transferred 
from the LS2 excitations to lower energies, and the ex- 
citation energy oi U + 2Jh is not observed in the ED 
results anymore. 

Here again the HF approximation reproduces all trends 
found in the ED results, but instead of the spectral weight 
transfer of the ED, one finds that the evolution of the HF 
wave functions simulates very well the ED findings at suf- 
ficiently large values of D. Systematic errors in the HF, 
due to the absence of quantum fluctuations, have been 
found only in the regime of small D < 0.1 eV, where the 
HF energies interpolate between the excitation energies 
U — Jh and U + Jh (lower by Jh than the corresponding 
ED values at 13 = 0) and the values expected from the 



FIG. 15. (Color online) Energy levels in the embedded bond 
(ij) along the c axis within the CG phase, see Fig. |4{b), as 
functions of increasing Coulomb interaction U for D — 30 
meV. Solid lines from HF and hollow circles from ED (the 
intensity scale and marker size indicate the spectral weight of 
the level). Other parameters as in set B in Table [D 



ED at sufficiently large D. Altogether, the agreement 
between the HF and ED is better in the regime of large 
D, where the quantum effects in the orbital space are 
suppressed. 



C. Localization by Coulomb interaction U 

So far, the results have shown that the HF approx- 
imation provides a realistic description of the electron 
density distribution in the occupied orbitals as well as of 
the excitation energies. This latter feature is not guar- 
anteed in the HF approach by itself, but follows here 
from the relatively large value oi U — 4 eV, which drives 
the electron localization and suppresses quantum fluctu- 
ations (except for spin fluctuations, fluctuactions related 
to the "pair hopping" term in Eq. (j2.16p . and the or- 
bital fluctuations along FM bonds). Here we show that 
the results obtained for the selected parameters are ac- 
tually valid in a broad range of values for the Coulomb 
parameter U. 

We fix Hund's exchange parameter Jh — 0.6 eV, which 
is responsible for the HS value {S = 1) at V'^+ ions in the 
ground state of the undoped system YVO3. Given this 
choice, we varied the parameter J7 in a quite broad range 
of values starting from U = 2.0 eVj2i As we have shown 
in Fig. [131 a relatively small value oi D — 0.03 eV suffices 
to suppress the orbital fluctuations along the considered c 
bond in the CG phase. In this regime of parameters, the 
excitation energies Em given by the HF approximation 
are in reasonable agreement with the ED ones (see Fig. 
[T5)) . and only the LS states have systematically lower 
energies than their ED counterparts, as discussed in Sec. 
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FIG. 16. (Color online) Occupations of spin-down pristine 
orbitals within the embedded bond {ij) along the c axis in the 
CG phase, see Fig. |3Ib), as functions of increasing Coulomb 
interaction U for D = 30 meV. Solid lines from HF and hollow 
circles from ED. Other parameters as in set B in Table H] 



As a matter of fact, this holds in the whole range 
oi U > 2.0 eV, i.e., in the entire physical range of local 
interactions relevant for transition metal oxides 

The dependences of the one-hole and the one-electron 
excitations on U are complemented by the electron den- 
sities projected onto the orbitals of t2g symmetry shown 
in Fig. [ini The c orbitals are almost filled in the entire 
explored range of U. They are relatively robust with re- 
spect to orbital polarization for the present small value 
of -D = 30 meV, but nevertheless the density gradu- 
ally decreases from riic = 1 at [/ = 2.0 cV to ~ 0.92 
at large U — 10 eV. This evolution of the electron densi- 
ties riic follows from the behavior of the orbital doublets 
{a, b} at both sites. 

First, in the regime of small U, electrons are not lo- 
calized as the charge excitation energy is the same as 
the hopping parameter. The ED gives then the densities 
0.58 (0.42) for the more (less) filled orbitals at the low- 
est considered value of U — 2.0 eV. The HF approach 
almost perfectly reproduces these densities. Increasing 
U suppresses the electron hopping along the bond, and 
the electron densities are gradually modified. The many- 
electron quantum state found in the ED is mainly de- 
termined by the fields acting on the orbital system: the 
JT effective fields from the neighboring atoms and the 
orbital polarization term. In the regime oi U < 3.5 eV 
the change of the electron densities is markedly slower 
in the ED than in the HF approximation, as here quan- 
tum fluctuations contribute. Only for higher values of 
U, both calculations agree with each other showing once 
more that the HF approximation is reliable in the entire 
strongly-correlated regime of large U due to the stabi- 



lization of symmetry breaking in spin-orbital space which 
suppresses quantum fluctuations. 



D. Ground state for the embedded bond 

One may wonder how accurate the HF approximation 
is for the total energy of the considered bond. When the 
bond forms, one gains binding energy due to the kinetic 
energy. The binding energy can be computed as the en- 
ergy difference between the energy of a bond and that of 
two single atoms. 



f(i7) = ^bo„d(C/)-2^atom(C/). 



(4.1) 



Thereby, one has to take into account the different num- 
bers of neighbors in each case. An atom is surrounded by 
six neighbors, each of them providing mean fields acting 
on the orbitals. In contrast, each atom of the bond has 
just five external neighbors, while the interactions along 
the bond are rigorously (approximately) included in the 
ED (HF). This generates a correction term a to Eq. (|4.ip 
included in Fig. 1171 where we analyze the superexchange 
contribution to the binding energy £ex{U) defined as 



£,ex(C/)=£({/)-£(^). 



(4.2) 



The binding energy £{U) (|4.ip calculated in the ED 
can be compared with the strong-coupling approach of 
Ref. |2^ to estimate the energy increments for the embed- 
ded bond and for the embedded atom in both magnetic 
phases, the CG and the GC phases. In the regime of 
large U ^ t, one finds asymptotic behaviors that suggest 
to use: 



U~3Jh 
9t 

U-\Jh 



(4.3) 
(4.4) 



as scaling parameters. The value of /3cg(C^) (|4.3p is ob- 
tained by considering HS excitations, while the value of 
Pgc{U) (|4.4I) is derived as an average value over the LS 
excitations allowed in this case. These quantities have 
been used to plot the superexchange energy £sex{U) (|4.2p 
for a single bond embedded in each magnetic phase con- 
sidered here, see Fig. [T71 We show the data for the same 
range of U values (2 eV < U < 10 eV) adopted in Figs. 
[T5]and[T6l respectively, which corresponds to /3(C/) < 1.0 
in both phases, and Pcg{U) > 0.024 (^Gc(t^) > 0.183) 
in the GG {GG) phase. 

The analysis of the electron density distribution in Sec. 
IIV CI suggests that the local Coulomb interactions are 
sufficiently strong to localize the electrons \iU > 3.5 eV, 
see Fig. [Tni and one expects that the strong-coupling 
expansion in powers oit/U should be valid. This is also 
confirmed by the present analysis of the binding energy 
£{U) (|4.ip — the numerical results obtained from the 
ED and the HF method agree very well in this regime 
for the GG phase, see Fig. 1171 Even more surprising 
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FIG. 17. (Color online) Hartree-Fock results (solid lines) for 
the superexchange contribution to the binding energy £sax{U) 
are compared to ED (circles) for an embedded bond (ij) along 
the c axis in the CG and the GC phases for D — 30 meV. 
For convenience the U dependence is expressed by a nonlinear 
scale (top) and by the perturbative scaling functions Pcg{U) 
(|4.3p and PcciU) (|4.4p . respectively. Other parameters as in 
set B in Table HI 



is the fair agreement between the ED and the HF cal- 
culations found for smaller values oi U < 3.5 eV, see 
Fig. [T71 This result shovifs that the strong orbital quan- 
tum fluctuations in the CG phase along the c axis^^ are 
well accounted for within the unrestricted HF approach. 
In the GC phase, the agreement is even better and the 
dependence on Pgc{U) is practically linear in the en- 
tire regime of considered values of [/, i.e., for U > 2.0 
eV. This confirms a quite strong electron localization in 
this phase, with the AF spin order preventing electronic 
transport and acting as a confining potential at each site 
of the bond. 



DEFECTS IN THE DILUTE LIMIT 



bors that form a cube shown in Fig. (TJ hence, the orbital 
degeneracy, which is controlled by the spin-orbital order 
in these compounds, is locally strongly affected by the 
defects4i The (occupied) orbitals in the neighborhood 
of the defect are polarized, i.e., rotated. This leads to a 
violation of the flavor conservation and, consequently, to 
a noticeable modification of the hopping matrix in the 
rotated t2g basis. Moreover, for each Ca-defect a hole 
appears in its proximity, bound to it by the Coulomb 
potential of the charged defect itself. Since these holes 
can move, controlled by the double exchange mechanism, 
they will affect the spin-orbital order further4^ In addi- 
tion, the chemical potential, due to the doped holes, will 
lie in a defect band inside the MH gap. Finally, also these 
defect states are subject to strong electronic correlations 
and display their own MH physics. 

Our main aim here is to investigate whether the MH 
multiplet splitting can be reliably obtained within the 
unrestricted HF framework for the 3- flavor t2g case, i.e., 
in presence of defect potentials, various orbital relaxation 
processes, kinetic energy of doped holes and electron- 
electron interactions. Here we consider the dilute limit 
of low electron doping x < 0.05, where the defects can 
be treated as well separated from one another. For more 
clarity, we also limit ourselves to hole doping for only one 
orbital flavor and one spin. Therefore, to demonstrate 
this, we focus below on several rather transparent cases 
of increasing complexity: 

(A) defect states in a multiband MH insulator in ab- 
sence of orbital polarization and of Coulomb inter- 
action among the electrons; 

(B) role of orbital polarization; 

(C) role of LR electron-electron interaction — Coulomb 
gap of defect states; 

(D) orbital polarization and relaxation of Coulomb gap 
of defect states. 

For the sake of clarity and simplicity, we will not con- 
sider here the effect of LR defect potentials, and the 
self-consistent screening of LR interactions, which results 
from the doping induced by the impurities. As our focus 
is on the strong-correlation aspects and in order to get 
easy-to-interpret eigenstates and spectral distributions, 
we shall not discuss the effects of disorder here, i.e., the 
defects will always be well annealed and in the dilute 
limit. 



A. Orbital polarization and the electronic structure 

The role of defects in transition metal oxides is sub- 
tle due to strong electron correlations, which manifest 
themselves in the multiplet structure of the MH bands. 
As a result of the orbital degeneracy, strong orbital po- 
larization and relevant relaxation processes occur in the 
vicinity of charged defects. Since a Ca-defect replaces an 
Y-ion in YVO3, each defect has eight vanadium neigh- 



B. Order parameter landscape and defects 

Unrestricted HF approach allows us to treat large inho- 
mogeneous systems, i.e., with defects, interfaces, nanos- 
tructures, etcetera. We have seen earlier for a two-flavor 
model that the Mott-Hubbard gap and the multiplet 
structure is faithfully reproduced for the strongly corre- 
lated transition metal oxide systems. A central require- 
ment of the unrestricted HF to work is broken spin and 
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FIG. 18. (Color online). Map of the orbital magnetization densities m^{k, I, m) (|5.1I) in the dilute well annealed CG phase with 
X = 0.015625 = 1/64, obtained for a 8 x 8 x 8 cluster (from positive to negative m^(k,l,m), see color scale on the far right). 
The coordinates A; = 0, . . . , 7 are along the horizontal axis, and m = 0, . . . , 7 along vertical axis, while the top row represents 
the / = and the bottom row the I = 1 plane. Four panels (from left to right) for different orbital flavor: 7 = c, 7 = 6, 7 = a, 
and 'y — a + b, respectively. Parameters: Vcc = 0, Vd = 1 eV and D — Q, others as in set B in Table H] 



orbital translational symmetry. Accordingly, vanadate 
perovskites, like Yi_a;Caa;V03, provide an ideal ground 
for this approach as spin and orbital orders in these com- 
pounds are present in a wide doping regime, actually up 
to the metal-insulator transition. Therefore, the unre- 
stricted HF is the method of choice in our case. Never- 
theless, it is not clear a priori how well the method works 
in the vicinity of defects in view of orbital polarization 
and hole motion. To explore this is our central goal. 

Moreover, the HF provides a basis for subsequent 
many-body perturbation treatments. However, it should 
be noted that here such a perturbative treatment will not 
be performed with respect to an 'uncorrelated' U = Q 
state, but with respect to the HF ground state. This 
latter state already features the MH gap and the higher 
multiple! structure. In fact, it was already pointed out 
that the fundamental MH splitting, i.e., the energy dif- 
ference of the LHB and the high-spin Hubbard band 
Ahs = {U ~ SJh) is obtained exactly within HF. We 
shall see below that there are other gaps, e.g. between 
defect states, which suffer from the well-known overesti- 
mation of gap energies within HF-theory. A subsequent 
many-body perturbation treatments of the screening of 
interactions is expected to cure such deficiences of the 
HF approach. 

In general, HF calculations require a self-consistent op- 
timization of all expectation values that define the final 
HF Hamiltonian, the HF ground state and the excita- 
tions. Complete information to construct the HF Hamil- 
tonian is thus contained in the order parameter land- 



scape, which is defined by the full set of effective fields. It 
contains already the rotations of occupied (and unoccu- 
pied) orbitals due to defects and doping. As an example, 
we display the diagonal occupation numbers in the CG- 
phase for the 3-flavor calculation performed with PBC in 
Fig. [TSl It shows maps of local magnetizations, 

m-y(k, I, m) — nj^{k, I, m) — n^^^k, I, m), (5.1) 

for an 8 X 8 X 8 cube of V ions and — from left to right 
— for the orbital character 7 — c,b,a and a + b. Two 
planes are presented: the top row represents the front 
face of the cube at / = 0, while the bottom row shows 
the / = 1-plane. 

The left panels reflect the C-type magnetic order of 
the electrons in c orbitals. The central panels show the 
same magnetic order for the b and a electrons. Here, in 
addition, one recognizes a checkerboard structure, which 
represents the G-type orbital order of the a and b elec- 
trons. In the right panel, magnetic densities of electrons 
in a and b type orbitals are summed up and reveal again 
the C-type magnetic structure. While the magnetization 
density has a perfect alternation between two sublattices 
in the / ~ 1 plane, four defects disturb it in the / — 
layer (top row), as seen for b and a orbitals, and also for 
their sum. 

The number of defects in the system corresponds in 
this case to hole doping x = 1/64. Here, the arrangement 
of the defects is well annealed, i.e., LR defect repulsion 
has been considered and, therefore, the defects are regu- 
larly spaced. The defects are visible in the top panel of 
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Fig. [18] and are located between the planes I = and 
1 = 1. The perpendicular feature in the a + 6-panel re- 
flects the motion of the doped holes along the c direction, 
but bound to the respective defects. From the 6-panel, 
one can also see that all 4 holes in the I = plane have 
6-down character, as assumed. We next analyze the one- 
particle excitation spectra for such a cluster obtained in 
the HF approach. 

C. Excited states in the 3-flavor model 

In the following, we shall investigate the HF excitation 
spectra for the four scenarios (A)-(D) introduced above. 
We will analyze only the CG phase as this phase is the 
relevant one for moderately doped Yi_2;Caj;V03, and we 
flrst focus on the dilute doping case oi x = 1/64. 

(A) We begin the discussion with the first scenario, 
where orbital polarization effects and LR Coulomb inter- 
actions are not considered, i.e., D = Q and Vcc = 0. The 
strong defect potential, Vd = 1.0 eV, confines a single 
doped hole per defect to the nearest neighbor V-sites of 
a defect, which form a cube. The corners of the Vs cube 
are equivalent with respect to the defect potential yet not 
with respect to the spin- and orbital- order. This implies 
that a {b J,) hole can delocalize only along a (vertical) 
bond along the c axis because of the AF order in the 
(horizontal) a6-plane. The kinetic energy of the a and b 
electrons in the a6-plane is quenched due to the double- 
exchange- type coupling to the spins of the c electrons, ^'^ 

To explore this symmetry breaking and its conse- 
quences for the excitations, we investigate in Fig. [19] 
the total DOS N{uj) together with the local partial DOS 
Nijai^^) foi' the two sites (0,0,0) and (0,0,1) of the ac- 
tive bond, i.e., the bond doped with a hole, and for a site 
(1,1,0) on an undoped spectator bond. The chosen spec- 
tator bond corresponds to the same spin polarization as 
the active bond, thus in principle the hole could delo- 
calize into its states too. However, this is not observed 
for the typical range of values of t (and other parame- 
ters). The total DOS in Fig. [lUJa), obtained for the 
8x8x8 cluster, shows the same multiplet splitting as 
discussed for the atom in Sec. IIIII We recognize the LHB 
and the three multiplet subbands forming the UHB. The 
fundamental MH gap opens between the LHB and the 
HS subband of the UHB; the splitting of these bands is 
given by the gap, 

Ahs = Ehs- Ei^HB = U -3Jh, (5.2) 

and is reproduced exactly in the HF approach. In con- 
trast, the actual MH gap, 

Amh = [/ - 3Jh - W,s, (5.3) 

is reduced by an effective bandwidth WeS- Due to the 
defects and their potentials, defect states are split off 
from the LHB. In other words, all single particle states 
in the vicinity of the defect are shifted upwards by the 



defect potential Vd (for simplicity, we consider here ony 
the case of a nearest neighbor defect potential). These 
states lie now inside of the MH gap. The transitions from 
the occupied defect states D to the upper HS Hubbard 
subband define the ingap absorption and the optical ab- 
sorption gap which we denote Aopt. At the same time, 
the spectral weight moves out of the LHB, similar as in 
the doped Hubbard modeL^^ 

The defect states would be completely filled, if there 
were not one doped hole per defect. This fixes the chem- 
ical potential fj, inside the defect band. It is useful to 
consider the integrated (averaged) electron density n{uj) 
per site, Eq. (|2.35|) . displayed in Fig. [T9][a), see right 
scale. The chemical potential fj, is then determined from 
the average of fi^ and /i+ which are obtained from the 
relation, 

n{fi^) = HQ — X ± e, (5-4) 

in the limit e — >■ 0, where no = 2 is the average number 
of electrons per site in the undoped 3-flavor model and x 
is the defect (hole) concentration. As can be seen from 
the inset in Fig. [TOT a). the chemical potential falls into a 
small gap at the upper edge of the defect band D. In the 
following, we denote this gap as the transport gap Atr as 
it would be a relevant measure for the conduction in the 
defect band. 

Deeper insight into the nature of the defect states and 
the transport gap is obtained by considering the local 
partial DOS N,^^{u}), shown in Figs. [I9];b)-[l9i;d) for 
sites (0,0,0) and (0,0,1), which form the active bond, and 
at the (1,1,0) V ion, a spectator site with two electrons. 
The defect, indicated by a red dot, lies in the center of the 
cube at (^, ^, ^). The two electrons at the spectator site 
(1,1,0) occupy the (c l) and {b l) local orbitals as can 
be seen from Fig. [TWd). They lie below the chemical 
potential /i indicated by the perpendicular dashed line. 

The doped hole has (6 |) character and is localized on 
the active bond displayed in Fig. ITW b) and (c) . On both 
sites the (c \.) orbitals are occupied, while the remaining 
(a \.) electron can delocalize on the vertical bond, i.e., 
parallel to the c axis, consistently with the FM correla- 
tions in the C-AF phase. This delocalization leads to 
a splitting into bonding and antibonding states on the 
active bond. 

Aba - 2^Vj\ + t^, (5.5) 

where the bonding-antibonding gap Aba is determined 
by the resonance integral t and the JT-interaction Vjt ~ 
Vab ~_\Vc in the case of the G-AO order. In Figs. ITWb) 
andlTOlc). Aba ~ 0.4 eV is essentially determined by the 
hopping matrix element t = 0.2 eV while Vjt ~ 0.02 eV. 
That the \a \.) occupation is asymmetric along the active 
bond is a result of the JT interactions and the underlying 
G-type orbital order. 

The origin of the 'transport gap' Atr is evident in this 
case and we identify it here with the gap between the 
topmost occupied orbital on a spectator bond, i.e., \b l) 
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FIG. 19. (Color online) (a) Total DOS N{uj) (left scale) and 
average electron filling n(aj) (right scale) in the dilute well 
annealed CG phase with defect concentration x = 1/64, ob- 
tained using an 8 X 8 X 8 cluster. The inset shows N{(jj) near 
the Fermi energy lj = 0. Panels (b)-(d) show individual DOSs 
Niacrioj) for orbital a and spin a, together with the total DOS 
N{uj) (thin lines). The atom coordinate, (0,0,0), (0,0,1), or 
(1,1,0), is indicated in each case and on the cube represented 
near each panel on the right. Parameters as in set B in Table 
m and 7? = 0, Voo = 0. 



FIG. 20. (Color onhne) (a) Total DOS N{ui) (left scale) and 
average electron filling n(aj) (right scale) in the dilute well 
annealed CG phase with defect concentration x — 1/64, ob- 
tained using an 8 X 8 X 8 cluster. The inset shows N{lj) near 
the Fermi energy uu = 0. Panels (b)-(d) show individual DOSs 
Niaaiy) for orbital ct and spin o, together with the total DOS 
N{uj) (thin lines). The atom coordinate, (0,0,0), (0,0,1), or 
(1,1,0), is indicated in each case and on the cube represented 
near each panel on the right. Parameters as in set B in Table 
B and D = 0.05 eV, Ke = 0. 
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orbital at site (1,1,0) and the unoccupied antibonding 
state on the active bond. This energy is 

^t.^ ^V^T + t^+ViT, (5.6) 

and for the present parameters one finds in Fig. [19] that 
Atr ~ 0.22 eV, and is essentially determined by the ki- 
netic energy. It is important to realize that within the 
HF-scheme occupied (unoccupied) orbitals correspond to 
electron removal (addition) energies. These processes 
need not to refer necessarily to the same defect, but can 
occur at different defects. Thus Atr is the energy gap for 
the hopping processes between different defects within 
the defect band. 

The energies of the occupied (c ],) orbitals differ 
strongly at different ions, see Figs. ITWbl-fTWd). This is 
a consequence of the different local occupation numbers, 
which enter the HF-energies. Obviously, this also affects 
the energies of the unoccupied orbitals. The defect states 
have their own MH gap as can be seen for example in Fig. 
fTOTd). Inserting an (a |) electron on site (1,1,0), in addi- 
tion to the already present (c |) and {h \) electrons, leads 
to the HS defect state |c | 6 | a |) at a; ^ 2.0 eV. Finally, 
at higher energies, one recognizes the LS defect multiplet 
states |c J, 6 J, 7 t), with 7 = a, &, c, respectively. 

(B) The effect of the Coulomb potential of the defect is 
to rotate the occupied orbitals in the vicinity of the defect 
such that their energy is minimized. This polarization ef- 
fect, described by iJpoi and controlled by the parameter 
_D, was not considered in the first scenario. The second 
scenario represents a calculation within the same param- 
eter set of the first scenario, with the exception for finite 
D = 0.05 eV. The total DOS iV(w) shows a clear change 
of the DOS of the defect states in the vicinity of /x, which 
is amplified in the inset of Fig. ISOT a') . While there was 
a pronounced transport gap in the absence of orbital po- 
larization, here this gap has almost disappeared (or is of 
marginal size). 

The main effects of orbital polarization can be revealed 
by inspecting the local partial DOS Ni^^{oj). As we ob- 
serve in Figs. [Sor b) and[20lc). which reflect the changes 
on the sites of the active bond, the (c' \.) orbital has now 
some admixtures of (a \.) and {h \.) flavors. More impor- 
tantly, Figs. I^CTb') and l^Uf c') show that the bonding and 
antibonding states are defined by linear combinations of 
I a 4,) and \b J,) states, and have no longer pure (a 4,) char- 
acter, as it was the case in Figs. [TWb) and ITW c). This 
mixing is a consequence of the fiavor non-conservation 
in the rotated t2g basis and of the resulting off-diagonal 
terms in the kinetic energy. 

The most important feature in this second (B) sce- 
nario is the formation of a bonding-antibonding splitting 
on the spectator bond, as can be seen in Fig. I20f d). At 
first glance, this may look puzzling as both states are oc- 
cupied, so there is no obvious net energy gain that may 
drive such a splitting. Inspection of energies and occu- 
pation numbers shows that the level-splitting is due to 
a combined mechanism: the rotation of the c orbital is 



favored when its energy coincides with the level of the 
bonding orbital. Thus, we observe here an orbital rota- 
tion induced by bonding-antibonding splitting on the spec- 
tator bonds. 

Interestingly, the antibonding band on the spectator 
bond is pinned to (but remains below) the chemical po- 
tential, whereas the antibonding state on the active bond 
is confined at (but remains above) the chemical potential. 
This leads to a vanishing transport gap Atr. However, 
this also implies a remarkable constraint in the system; 
namely, the hole cannot delocalize from the active bond 
into the opposite spectator bond with the same spin ori- 
entation. 

Another quite important consequence of the orbital fla- 
vor mixing is the collapse of the charge polarization on 
the active bond. Whereas in the absence of orbital polar- 
ization in scenario (A) there is a pronounced accumula- 
tion of a-type charge at site (0,0,1) and of hole-density at 
site (0,0,0) triggered by the JT interaction, one observes 
an essentially complete collapse of the combined a plus 
b charge polarization on the active bond in the case of 
D = 0.05 eV, see Figs. M^) anddOtc). 



D. Role of long-range Coulomb interaction 

Next we turn to the effects of the LR Coulomb in- 
teraction in presence of defect states. In insulators the 
Coulomb interaction is screened by the dielectric con- 
stant eo of the system; yet, the interaction keeps its 
^ l/eor LR character. Here Eq represents the back- 
ground dielectric screening due to the 'core'-electrons. 
These are in our case all electrons apart from those of 
^2g-type. The screening arising from t2g electrons, which 
is for instance essential for the MIT, is explicitly included 
in the Hamiltonian of the system. 

In the next scenario (C), the LR Coulomb interaction 
Eq. (j2.20l) is taken into account and parametrized via 
Vcc (instead of Eq), while the orbital polarization is ne- 
glected {D = 0), as in scenario (A). The DOS is displayed 
in Figs. l^lTa) andl^TTb'). Its comparison with Fig. [TWa) 
shows that the LR Coulomb interaction has no signifi- 
cant effect on the ionic multiplet structure, as one might 
expect. There is, however, a significant change in the 
size of the energetic splittings between the different de- 
fect states in the vicinity of the chemical potential, as 
one can see from the insets in Figs. [T9l a) and [2lT a'). 
The fine structure of the levels is here basically as in 
Figs. [TWc)-fTWd). The peaks in the inset stem from 
the occupied {(c 1), (6 \.)} orbitals at the spectator sites 
as well as from the (a I) bonding/antibonding levels at 
the active bond. The main effects of the nonlocal LR 
electron-electron Coulomb interaction are: (i) to enlarge 
the bonding-antibonding splitting, and (ii) to increase 
the polarity of the active bond. 

The transport gap Atr opens between the topmost oc- 
cupied |c 4 & 4-) states on the spectator site and the unoc- 
cupied I a I) antibonding state on the active bond. Fig- 
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FIG. 21. (Color online) (a) Total DOS iV(aj) (left scale) and 
average electron filling n[uj) (right scale) in the dilute well 
annealed CG phase at hole doping x = 1/64, obtained using 
an 8 X 8 X 8 cluster. The inset shows N{ijj) near the Fermi en- 
ergy uj = Q. Calculation includes the LR Coulomb interaction 
Kc = 0.2 eV. (b) Dependence of DOS N{ijj) of defect states 
and Coulomb gap near the Fermi energy oj = on increasing 
strength of LR Coulomb interaction: Vco = (solid, black), 
0.2 eV (dashed, red), 0.5 eV (dashed-dotted, blue), 1.0 eV 
(dashed-double-dotted, dark green). For increasing strength 
of 14e the features in N{uj) move away from the Fermi energy 
u) — 0. Other parameters as in set B in Table I and D — 0. 



ure HHb) shows that it increases almost linearly with 
increasing LR electron-electron interaction < Vee < 1.0 
eV. Indeed, this dependence of Atr on Vcc is summarized 
in Fig. [2H The line A-B shows the evolution of Atr 
with the increasing value of Vcc for vanishing polariza- 
tion field D ~ 0. In the absence of nonlocal electron- 
electron Coulomb interaction, the origin of the gap is the 
bonding-antibonding splitting due to the kinetic energy 
gain of the single a-electron on the active bond. With 
increasing Vcc, the charge polarization of the bond grows 
and the kinetic energy is suppressed. At large Vee, the 
gap Atr is determined by the nonlocal Coulomb interac- 
tions, i.e., essentially by the nearest- and further neighbor 




FIG. 22. (Color online). Transport gap Atr in the well an- 
nealed CG- AO phase with doping concentration a; = 1/64 as 
a function of electron-electron interaction coupling Ko and or- 
bital polarization strength D. Results obtained for an 8 x 8 x 8 
cluster with the parameters as in set B of Table I. 

ones. The increase of Atr along the line A-B in Fig. [22] 
may be approximated in terms of the electron-electron 
interaction parameter Vcc, the hopping integral t and the 
JT coupling VjT as follows, 

Atr « ViVec/^ + Vjt)^ + i2 + Vcc/2 + VjT ■ (5.7) 

We emphasize that, in these systems, the Coulomb 
gap in the defect states is a consequence of the com- 
plex structure of the defects in combination with the 
electron-electron interactions. We note that this mecha- 
nism is distinct from the disorder-induced Coulomb gap 
of Shklovskii and Efros.^^ We have already seen that the 
doped holes in the dilute doping regime have no effect 
on the multiplet splitting away from the defects. On the 
other hand, on general grounds we would expect that the 
present Coulomb gap in the defect band is overestimated 
within the HF approximation, and will be reduced by 
dielectric screening, i.e., obtained by a many-body treat- 
ment beyond the HF approach. We shall see next that 
already the inclusion of orbital polarization gives rise to a 
substantial screening and reduction of the Coulomb gap. 

Finally, we consider scenario (D) that includes both LR 
electron-electron Coulomb interaction and orbital polar- 
ization effects in the vicinity of the defects. The DOS is 
displayed in Figs. I^HTa) andl^Hlb). Comparison of the in- 
sets in Figs. BHT a) and l^lT a) shows the strong reduction 
of Atr in the present case oi D — 0.05 eV from the value 
found with D = 0. Similarly, Figs. [M^b) andHHb) make 
apparent the strong reduction of the gap at D — 0.05 eV 
not only for Vcc — 0, but also for larger values of Vcc- 

To summarize, the central observations concerning the 
transport gap are: (i) the decrease of Atr with increasing 
polarization D, and (ii) the increase of Atr with increas- 
ing LR electron-electron Coulomb interaction strength 
Vcc- These findings are summarized in Fig. 1221 The 
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FIG. 23. (Color online) (a) Total DOS N{uj) (left scale) and 
average electron filling n(aj) (right scale) in the dilute well 
annealed CG phase with x = f/64, as obtained using an 8 x 
8x8 cluster for finite LR Coulomb interaction 14c = 0.2 eV. 
The inset shows N{lj) near the Fermi energy uu = 0. (b) DOS 
N{uj) of defect states and Coulomb gap near the Fermi energy 
uj — for increasing strength of LR Coulomb interaction: 
Voo = (solid, black), 0.2 eV (dashed, red), 0.5 eV (dashed- 
dotted, blue), 1.0 eV (dashed-double-dotted, dark green). For 
increasing strength of Vcc the features in N(lu) move away 
from the Fermi energy lj — 0. Other parameters as in set B 
in Table I and D = 0.05 eV. 



strong impact of the orbital polarization D-term on Atr 
is seen between the A and C points, where the transport 
gap is reduced by a factor ~ 3 when D increases from 
to 0.05 eV. This strong suppression of the transport gap 
change found here for relatively weak polarization inter- 
action D is consistent with fast orbital rotation induced 
by this term, presented in Sec. IIIIl 



VI. DISCUSSION 
A. Gaps in the electronic structure 

In our HF study we have encountered three fundamen- 
tal energy scales typical of doped MH insulators. They 
are displayed in Fig. [23] for the dilute doping regime 
X < 0.05 investigated here. In decreasing-size order, 
these energies represent: (i) the Mott gap computed as 
the energy of the HS charge transitions Eq. (|5.2p . (ii) 
the optical gap, Aopt, relevant for the in-gap absorption, 
which opens between the occupied (unoccupied) defect 
states inside the MH gap and the unoccupied (occupied) 
states of the HS part of the UHB (the LHB), and (iii) the 
lowest energy excitation gap within the defect states Atr 
Eq. (|5.7p . of relevance for transport. The reported data 
are obtained for a finite value of the LR Coulomb inter- 
action strength Vcc — 0.2 eV, and for two different values 
of the orbital polarization parameter D, corresponding 
to scenarios (C) and (D). Note that the gaps Ahs and 
Aopt in Fig. [M] are inferred from the peak positions in 
the DOS and correspond to the absorption maxima and 
not to the onset of the absorption. 

It is evident that only Atr is significantly affected by 
the orbital polarization term (Fig. B^ . We note here 
that while the screening of electronic interactions is gen- 
erally obtained by many-body perturbation theory be- 
yond the HF approximation, for example by resumma- 
tion in random phase approximation, the screening of 
Coulomb interactions at defects is already, to a large ex- 
tent, contained on the Hartree level and taken care of by 
charge and orbital relaxation effects around the defects. 
This explains why finite orbital polarization D leads to 
such a strong reduction of the transport gap. 

These characteristic scales may be compared with 
the experimental data^ for the doping dependence of 
the optical absorption spectra for Yi_a;Caa;V03 and 
Lai_a;Sra;V03. They reveal HS Mott transition energies 
Ahs of about 2.2 eV and 1.8 eV, respectively. The Mott- 
gap energy is essentially independent of doping, as found 
in the present HF approach. Only close to the MIT, the 
intensity of the MH absorption fades away. The most 
remarkable feature in optical spectroscopy is the appear- 
ance of an in-gap absorption related to defect states: the 
intensity of this absorption grows with doping. The char- 
acteristic energy of the in-gap absorption Aopt (read off 
at the peak positions) is independent of doping in the 
dilute limit. 

At higher defect concentrations, which go beyond the 
present model calculations designed for the dilute limit, 
i.e., for X < 0.05, the energy of the absorption peak 
decreases and eventually collapses for doping concen- 
tration approaching the MIT. For example, in the case 
of Yi_^Ca^V03, Aopt ~ 1.2 eV up to x = 0.10.^ In 
Lai_a;Sra;V03, this absorption is centered at about 0.80 
eV for X < 0.10, while the in-gap absorption peak shifts 
to ~ 0.3 eV at the doping concentration x = 0.168, i.e., 
close to the MIT4 In the high concentration regime it 
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FIG. 24. (Color online) Three characteristic gaps in the spec- 
tra for increasing doping x: the Mott gap Ahs Eq. (|5.2p 
deduced from the energy of the HS transitions (triangles), 
the optical gap Aopt between the defect states and the HS 
UHB (squares), and the transport gap Atr within the de- 
fect states (circles), determined from the flat region in n(tj) 
at a; = 0. The three gaps are obtained for the dilute well- 
annealed doped CG phase using an 8 x 8 x 8 cluster, and for 
two values of D: D = (solid lines, black) and D = 0.05 eV 
(dashed lines, red). Parameters: as in set B in Table I and 
Voc = 0.2 eV. 



will be important to include the dielectric screening due 
to the doped holes. 

Activated transport in Lai^j^SrajVOs has been re- 
ported first by Dougier and HagenmuUeij^who observed 
an activation energy Atr ^ 80 meV at x ~ 0.05 dop- 
ing. Transport in these systems has been extensively dis- 
cussed by Motl3 and identified as Anderson type, i.e., 
controlled by defects. For Yi-a^Caa^VOa, activated be- 
havior of the resistivity was reported for the Ca-doping 
range x < O.SJ^i^^ In particular. Sage et al^ reported 
the activation energies Ea = 0.124, 0.106, and 0.064 eV 
at doping x = 0.1, 0.2, and 0.3, respectively. The val- 
ues for Atr in Fig. [53] are of similar size when orbital 
polarization D = 0.05 eV is included. Noticing these 
qualitative trends, we point out here that a quantitative 
analysis of the above experimental data will be possible 
only after extending the present model to the regime of 
higher doping x > 0.1. 



B. Orbital density distribution 

Spin-orbital order in RVO3 compounds is typically 
discussed by considering only the subspace of a and b 
orbitalsJ^iS While such a simplified picture may be suf- 
ficient for the undoped compounds, it becomes question- 



FIG. 25. (Color online) Average total electron filling n-, in the 
c orbitals (ric, circles) and in the {a, b} orbital doublet {uab, 
squares), as obtained using an 8 x 8 x 8 cluster in the dilute 
well annealed CG phase, for increasing doping x and for two 
values of the orbital polarization interaction: D = (solid 
lines), and D — 0.04 eV (dashed lines). Other parameters as 
in set B of Table I, and Ke =0. 



able when the system is doped. One important and ob- 
vious perturbation mixing the orbitals is the interaction 
with the dopands. In fact, this interaction triggers the or- 
bital polarization i/poi, induced by the polarization con- 
stant D in Eq. (|2.12l) . Figure [25] displays the total filling 
within c and a + b orbitals (defined with respect to the 
original unrotated orbital basis) as function of doping 
X. We see that, as expected, the c occupation remains 
unchanged when the orbital polarization is absent (at 
D = 0), i.e., holes go only into the a and b orbitals. In- 
stead, the situation is reversed for D = 0.04 eV and the c 
occupation is here more strongly reduced than the a + b 
occupation. 

This, at a first glance, very surprising and counter- 
intuitive result can be understood by recalling that the 
lowest (occupied) local states rotate from pure c char- 
acter at Z) = to a linear combination of c and {a, b} 
at finite D in the vicinity of the defects, see Sec. IIIII 
It is important to realize that this rotation involves all 
V-neighbors of the defect. Thus, the concentration of V- 
ions affected by the redistribution of electronic density 
is much larger than the concentration of doped holes — 
this effect explains the strong change of occupancy seen 
in Fig. [15] 

Considering the strong admixture of a and b orbital 
character in the rotated c'-orbital, one may wonder 
whether this does not necessarily imply a complete break- 
down of any analysis based on a two-flavor description. 
Such concerns are certainly justified when matrix ele- 
ments of local orbitals come into play as, for example, in 
the calculation of the intensities of the transitions inves- 
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FIG. 26. (Color online) Occupied part of the (spin-orbital)- 
resolved defect state DOSs Niacrii-o) for orbital a and spin a at 
the atom (1,1,0) which belongs to undoped {spectator) bond, 
as in Fig. [20jd), obtained in the dilute well annealed doped 
CG phase {x = 1/64) using an 8 x 8 x 8 cluster. Parameters 
as in set B of Table |fl D = 0.05 eV and Kc = 0. 



tigated in optical or RIXS spectroscopies. We emphasize 
that in cases where models merely rely on the fact that 
the rotated c'-orbital is the lowest orbital in the t2g sec- 
tor and is always occupied, qualitative conclusions based 
on two-flavor models remain valid. 

Further insight into the variation of electronic filling in 
different orbitals at increasing doping x can be gained by 
inspecting the occupied single particle states at a specta- 
tor site (1,1,0) in the vicinity of a defect at (i, i, i) for 
finite orbital polarization D = 0.05 eV, as displayed in 
Fig. [26l This figure highlights the formation of bonding 
and antibonding states on spectator bonds, which involve 
the (a I) and (6 \) orbitals as well as the (c ],) orbital 
(but with a much reduced splitting). An important fea- 
ture here is that the chemical potential is pinned at the 
upper edge of the antibonding band. This guarantees 
that the spectator states are blocked for doped holes! 



Another important effect of D is that, with respect 
to the rotated basis, the hopping matrix no longer con- 
serves flavors, which is the case for the original ba- 
sis. Therefore, large polarization D leads to a strong 
mixing of orbitals with different character in the vicinity 
of defects also via the kinetic energy. One may won- 
der which mechanism leads to the bonding-antibonding 
splitting on the spectator bond although both states are 
occupied such that there is no evident energy gain as, for 
instance, in the case of the JT splitting or the Peierls dis- 
tortion. Here, we observe that the splitting is induced by 
the orbital rotation (at finite D) and the appearance of 
the off-diagonal hopping processes along the (vertical) c 
axis. Remarkably, even for the rotated c orbital, hopping 
along the c axis is now possible, as one can see from a 
small splitting of the DOS for c electrons in Fig. [551 



VII. SUMMARY AND CONCLUSIONS 

Before summarizing, we recall the fundamental prob- 
lem in strongly correlated systems: the typical multiplet 
splitting of the transition metal ions can be obtained in 
the localized limit by a Hartree factorization performed 
with respect to an optimal local basis set, i.e., a basis with 
occupation numbers being close to either 1 or 0. That 
is, the local spin-orbital states should be either occu- 
pied or empty — intermediate values for the occupation 
numbers would signal the breakdown of the Mott gap 
in the Hartree approximation, which ultimately happens 
when the kinetic energy increases and dominates over 
electron localization triggered by the local interactions. 
The former requirement, however, leads to complications 
in the vicinity of defects. This follows from a rotation 
of the occupied orbital states of the V-ions due to the 
local orbital-polarization fields associated with charged 
defects. 

In principle, a local rotation could be specified that 
removes the off-diagonal polarization terms and defines 
a new optimal local basis, i.e., again with occupancies 
either close to 1 or 0, and a fully developed Mott gap. 
This simple scheme is not applicable when some nonlocal 
terms are present (as e.g. the kinetic energy term). Yet, 
it is important to recognize that the unrestricted Hartree- 
Fock (HF) (and not Hartree) scheme, i.e., including the 
relevant off-diagonal contributions from (local and non- 
local) interactions, can be used to determine the Mott- 
Hubbard split bands in the presence of defects. We have 
found that this also holds true for the three-fiavor case. 
In particular, we have shown, by comparison with exact 
diagonalization, that the HF approach provides a sur- 
prisingly faithful description not only of the occupied, 
but also of the unoccupied higher multiplet states. It is 
the latter aspect that is particularly surprising here as it 
is known that the HF approach is designed to describe 
the occupied states in the best possible way, while for 
the unoccupied states one expects a considerably poorer 
description (e.g. electron affinities). 

Basic features of the character of wave functions and 
excitation spectra, which reflect the interplay of strong 
correlations and the orbital polarization fleld D due to 
a defect, are already well represented by the HF study 
of a single ion or a bond embedded into the electronic 
structure of a large cluster reflecting the Mott insulator 
with spin-orbital order. These studies and their compar- 
ison with exact diagonalization allow us to conclude that 
the HF calculation describes well: (i) the effect of the 
orbital polarization on the wave functions in the three- 
flavor case, and (ii) the multiplet structure for an atom 
(Sec. IIII[) . A further investigation for an embedded bond 
(Sec. IIVP highlights the interplay of the orbital polar- 
ization field D and the kinetic energy. Comparison of 
the HF results for partial densities and excitation spec- 
tra (i.e., photoemission and inverse photoemission) with 
exact diagonalization data shows good agreement in the 
whole D range. An exception are extremely small values 
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of D, where the exact wave function is entangled, a fea- 
ture that cannot be captured within the HF. A particu- 
larly good agreement between the HF and exact diagonal- 
ization is obtained for the binding energy basically for all 
values of J7 > 2 eV. The binding energy for the CG phase 
is found to be much larger than that for the GC phase. 
This higher binding energy, together with the double ex- 
change mechanism realized on ferromagnetic bonds, is 
responsible for the increasing stability of the GG phase 
on increasing the doping and leads to the experimentally 
observed magnetic phase transition in Y1_2.Ca3.VO3 at 
very low doping, x ^ 0.02.^» 

The description of doped carriers in Mott insulators, 
such as in the Ri-^^x^O^ compounds {R—Y, La and 
A^Ca., Sr, etc.) or in LaVOs/SrVOa superlattices^ 
is a computational challenge due to the interplay of 
strong correlations, defect-induced local deformations of 
the wave functions and spin-orbital order. It is evident 
that to deal with strong correlations and defects, for suf- 
ficiently large systems, requires carefully chosen approx- 
imations. In this work, we have shown that the unre- 
stricted HF method fulfils all essential requirements. In 
particular, we have seen that even the unoccupied states 
at high energy, i.e., those reflecting the multiplet struc- 
ture of the Mott insulator at orbital degeneracy, are well 
described in a doped system. The same holds true for 
the perturbed electronic structure in the vicinity of the 
charged defects. The HF approach is an efficient scheme 
which maps the interacting electron problem onto the 
problem of a single particle moving in a self-consistently 
determined field of the other electrons. Therefore, it is 
applicable for large systems, systems with defects, het- 
erostructures and interfaces. Very relevant for the suc- 
cess of the method in describing Mott insulators are the 
broken symmetries reflecting the underlying spin-orbital 
order. Fortunately, the spin-orbital order in Ri-xAxYO^ 
compounds persists up to the metal-insulator transition, 
thus the assumption of spin-orbital order in the unre- 
stricted HF scheme is here fully justifled. 

Having shown that the HF approach provides a faith- 
ful description of excited states of single defects in sys- 
tems with broken symetries, we investigated the elec- 
tronic structure of an 8 x 8 x 8 cluster for flnite defect 
concentrations in the dilute doping regime. The aim here 
is to resolve the complex nature of defect states and their 
manifestation in the excitation spectra. A single A de- 
fect, for example a divalent A ion in Ri-^A^YOa, adds a 
hole into the t2g shell of V(3d) states. In the dilute limit, 
this hole is not freely moving but is pinned by the charge 
of the defect essentially to a cube of V ions surrounding 
the defect A in its center. 

The present study demonstrates that the excitation 
spectra of the intrinsic, orbital-degenerate Mott insula- 
tor are strongly modified by adding defects. We address 
here only some points relevant for experimental studies. 
An important observation is that the fundamental exci- 
tations, such as the Mott-gap and the multiplet energies, 
are not affected by doping in the dilute regime. The 



defects generate defect states inside the Mott-Hubbard 
gap, as observedi^ Due to the intrinsic doping, the chem- 
ical potential lies inside the defect band. The spectral 
weight of defect states is taken from the LHB. We have 
argued that the in-gap absorption observed in the opti- 
cal conductivity experiments performed for several doped 
vanadate systems can be identified with these states. 

Moreover, we observe in the HF excitation spectra for 
the GG phase a transport gap inside the defect band, 
which corresponds to the removal of an electron at one 
defect and the addition at another defect in the neigh- 
borhood. In the absence of electron-electron interac- 
tions, this gap is essentially determined by a bonding- 
antibonding splitting resulting from the c-axis kinetic en- 
ergy of a doped hole confined by the defect. The inclu- 
sion of LR Coulomb interactions between electrons ^ Vcc 
leads to an approximately linear increase of the trans- 
port gap with Vac- Thus, these interactions promote the 
transport gap to a Coulomb gap. Combining Coulomb 
interactions and the mechanism of orbital polarization we 
find a reduction of the size of the Coulomb gap. This can 
be considered as screening in the vicinity of the localized 
defect, i.e., screening contained already at the HF level. 

An important motivation of this work was the study 
of the three-flavor case, particularly in connection with 
defects in the RVO3 perovskites. Usually spin-orbital or- 
der in the perovskite vanadates is discussed in terms of 
a simplified two-flavor ({a, b}) model4^ We have shown 
here that defects lead indeed to a strong change of occu- 
pied orbitals, e.g. the occupied c-orbitals, in the vicin- 
ity of charged defects due to orbital polarization. Yet, 
we have also shown that the topmost occupied local or- 
bitals, i.e., the orbitals that are relevant for the doped 
holes, are mainly of {a, 5}-character as in the two-flavor 
description. It is this latter observation that leads to the 
conclusion that the interaction of doped holes and the 
spin-orbital degrees of freedom are in fact similar in the 
two models. 

Summarizing, the central result of this paper is es- 
tablishing that the unrestricted HF approach is well de- 
signed to describe the electronic structure of charged de- 
fects in doped transition metal oxides with active orbital 
degrees of freedom, in presence of strong electron correla- 
tions. The present study provides valuable insights into 
the changes of the electronic structure in doped Mott- 
Hubbard systems under doping. While the positions of 
the Hubbard subbands are not affected by doping, their 
spectral intensity changes and new defect states occur 
within the Mott-Hubbard gap. These new states are ob- 
served in the optical spectroscopy and their weight in- 
creases with doping. 

Finally, we remark that the HF approach represents a 
natural basis for a subsequent many-body perturbative 
treatment. It should be noted that the unrestricted HF 
method already provides the Mott-Hubbard gap and the 
higher energy multiplet excitations. Thus, the proper 
many-body treatment would start from a reasonable 
Slater determinant, which reflects the correlated elec- 
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tronic structure of the spin-orbital ordered Mott insu- 
lator. We have shown that the fundamental Mott gap 
is correctly described by the HF in the framework of 
the multiband Hubbard model. Therefore, the many- 
body treatment will hardly affect the Mott-Hubbard gap, 
yet it may further improve higher multiplets, and it will 
certainly contribute to the screening of the long-range 
Coulomb interactions and lead to an extra reduction of 
the Coulomb gap in the defect states. 
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Appendix: Hartree versus Hartree-Fock 
approximation and optimal basis 

The remarkable successes of the LDA+U approach^^ 
in describing the various aspects of the electronic struc- 
ture of strongly correlated materials ^^"-^ including 
the magnetic order and the MH gaps, suggest that the 
Hartree (or mean field) approximation for the electron- 
electron interactions in the appropriate orbital basis 
could be sufficient to design an efficient and realistic 
scheme to determine the electronic structure. This ex- 
perience, however, is based on undoped systems, where 
the occupation numbers for the atomic orbitals (such as 
t2g or Eg orbitals at transition metal ions in the corre- 
lated oxides) used in the LDA+U calculations are either 
close to 1 or close to 0. This latter condition is actu- 
ally necessary to produce an essentially correct multiplet 
splitting, i.e., these orbitals are already the properly cho- 
sen orbitals to implement the electron interactions in the 
Hartree scheme. This is always the case in an undoped 
Mott insulator, which explains the success of the LDA+U 
approachj^^^-iS^ Note also that the orbital basis is then 
identical at every site. 

However, in a system with defects the occupied orbitals 
(here we consider J,-spin ones) belong to an orbital basis, 
which adjusts itself due to the external polarization field 
(|2.12[) . see Sec. IHII As a result, for an atom one finds new 
rotated orbitals {|Cn)}, introduced in Sec. IHII Similar 
situation occurs in a bond but in addition the orbitals 
are then delocalized over two atoms in the CG phase, 
see Sec. IIVI Hence we write the HF orbitals for ^-spin 
electrons {|Cn)} as a linear combination of the original 
t2g orbitals, {\ia), \ib), \ic)} (|2.6p at site i, 

i 7— a 



TABLE IL Decomposition of the local orbital basis Eq. 
HA. II) onto the original t2g basis given by the elements lai""*!^ 
(|A.2|I for ^-spin electrons, as obtained for an atom (atom, 
n = 1,3,5) and for a bond with two atoms i = 1,2 (bond, 
n = 1, ■•• ,6) in the CG phase, with D = 0.05 eV. Other 
parameters as in Table I: set A (atom) and set B (bond). 







atom 




bond 




orbital 
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i = 2 


total 
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0.8616 


0.0058 


0.0058 
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0.0832 
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0.4942 




a 


0.0552 
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0.4942 


^2 


c 


— 
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0.0430 


0.0232 
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0.3293 


0.1611 


0.4904 
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0.3293 


0.4904 
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0.0786 
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0.3063 


0.4213 
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0.4861 


0.3063 


0.1150 


0.4213 
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c 




0.0605 


0.0605 


0.1210 




b 




0.1470 


0.2925 


0.4395 




a 




0.2925 


0.1470 


0.4395 



Here i = 1 and n = 1,3,5 for an atom, while i = 1, 2 and 
n — 1, . . . , 6 for a bond along the c axis. As an example, 
we give in Table [III the expansion coefficients, 

|ag)p^|(e„N7)P (A.2) 

obtained for a representative value of 13 = 0.05 eV. These 
coefficients were obtained by applying the unrestricted 
HF approximation as described in Sees. IIHI and ITVl 

It is important to realize that the Fock terms that stem 
from interorbital Coulomb interactions between electrons 
with the same spin play an essential role in optimizing the 
orbital basis near charge defects. Indeed, the polarization 
term (|2.12p acts as a field and contains on-site interorbital 
averages, {cl^^c^^^} with a /3, which are responsible for 
the orbital rotation and the change of the local basis. By 
comparing the occupied orbitals for an atom and for a 
bond in the C-AF phase at finite value of i = 0.2 eV, see 
Table [Hi one finds that further change of the orthogonal 
basis occurs at each site due to finite kinetic energy along 
the FM bond Note that the orbital labels {n} in Table 
|TT]are assigned to the orbital functions {^n} in sequence 
of their increasing HF energy. The orbitals |^i) and j^s) 
are occupied for an atom, while for the bond the number 
of occupied orbitals is doubled and also the orbitals |^2) 
and 1^4) are occupied — they are characterized by similar 
orbital densities as those in the orbitals |^i) and j^a) for 
the atom, respectively, but the energies of the latter two 
orbitals are interchanged in the bond. 

The local electron-electron interactions are rotation- 
ally invariant, and the Hartree scheme could in prin- 
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ciple be applied to any set of locally orthogonal orbitals 
including the basis of rotated orbitals discussed above. 
If one could find such a rotated basis with the self- 
consistently determined orbital occupations being just ei- 
ther or 1, Hartree approximation would again be again 
sufficient as in the LDA+J7 scheme. Mean fields would 
then suffice to come quite close to the exact solution (as 
for a single atom, see Sec. IIIip . and to reproduce the mul- 
tiplet structure of the UHB. In such a case, Fock terms 
simply vanish and one can neglect them from the very 
beginning despite the fact that terms that could drive fi- 
nite values of the off-diagonal averages are present in the 
HF Hamiltonian. 

Summarizing, we conclude that the off-diagonal Fock 



elements are essential within the HF scheme and are nec- 
essary to arrive at the final optimal orbitals adjusted to 
the defect states. It is only in this orbital basis that the 
electronic structure in a Mott insulator faithfully repro- 
duces the multiplet splittings. Therefore, possible exten- 
sions of the LDA+U scheme to doped Mott insulators 
will have to use downfolding procedure^*^^ to a proper 
tight-binding model^°^ in which full unrestricted HF cal- 
culations could be performed. Designing such a scheme 
which has to go beyond the present LDA+U approach^ 
is indeed very challenging — it would make it possible 
to investigate not only doped Mott insulators but also 
interfacesr^i^ heterostructures^iSS, or other composite 
materials with defect states in correlated insulators in 
the future. 
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